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Abstract 

This  report  presents  a  summary  of  the  research  supported  by  the  AFOSR  during  the  period 
February  1  2005  to  May  31  2008  on  the  formulation  and  analysis  of  integration  algorithms  for  the 
nonlinear  dynamics  of  solids  and  structures.  The  project’s  main  focus  has  been  the  development  of 
energy-dissipative  momentum-conserving  time-stepping  algorithms  (or,  in  short,  EDMC  schemes) 
for  finite  strain  plasticity  and  nonlinear  coupled  thermoelasticity.  These  schemes  lead  to  numerical 
solutions  that  exhibit  exactly,  by  design,  the  conservation  laws  of  linear  and  angular  momenta,  as 
well  as  the  exact  physical  dissipation  characteristic  of  these  inelastic  systems.  This  latter  property 
leads  to  a  much  improved  numerical  stability  when  compared  with  classical  schemes,  avoiding  in 
particular  their  observed  numerical  instabilities  in  the  nonlinear  range.  The  design  of  the  new 
algorithms  is  based  on  a  complete  mathematical  analysis  of  the  discrete  dynamical  system,  with 
the  above  conservation/dissipation  properties  incorporated  and  proven  rigorously  and  for  general 
models  and  conditions  (e.g.  independent  of  the  time  step).  In  addition,  the  observation  of  the 
crucial  role  also  played  by  the  spatial  discretization  has  led  to  the  development  of  new  locking-free 
assumed  strain  finite  element  methods  for  the  implementation  of  energy- moment um  schemes. 


KEYWORDS:  nonlinear  dynamics,  energy-dissipative  moment  inn- conserving 
algorithms,  finite  strain  plasticity,  nonlinear  coupled  ther¬ 
moelasticity,  assumed  strain  finite  element  methods. 
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1.  Objectives 

The  research  considered  in  this  grant  addressed  different  aspects  of  the  numerical  analysis 
of  the  dynamics  of  nonlinear  solids  and  structures.  The  main  goal  was  the  development 
of  stable  time-stepping  algorithms  for  nonlinear  dynamics.  The  focus  was  on  inelastic 
solids,  including  finite  strain  elastoplastic  and  coupled  nonlinear  thermomelastic  models 
of  three-dimensional  solids  and  of  structures  (like  shells).  In  addition,  the  project  has 
considered  the  first  steps  in  the  extension  of  these  ideas  to  the  modeling  of  failure  in  solids 
and  structures  in  the  dynamic  range,  through  the  consideration  of  finite  elements  that 
incorporate  strong  discontinuities  for  the  modeling  of  cracks  and  other  similar  localized 
solutions. 

This  effort  significantly  extended  the  range  of  application  of  the  methods  that  we  have 
developed  previously  with  the  support  of  the  AFOSR  for  elastic  problems.  In  this  way, 
the  main  goal  can  be  summarized  as  the  development  of  time-stepping  algorithms  that 
rigorously  exhibit  the  non-negative  energy  dissipation  characteristic  of  these  inelastic  sys¬ 
tems,  while  preserving  the  conservation  laws  of  linear  and  angular  momentum  and  the 
associated  relative  equilibria  of  the  dynamical  system.  Existing  classical  techniques  were 
developed  in  the  context  of  linear  problems  and  do  not  exhibit  these  crucial  properties 
when  employed  in  nonlinear  applications.  In  fact,  classical  time-stepping  algorithms  (like 
New  mark  schemes  and  their  variations)  have  been  observed  to  numerical  instabilities  in 
the  nonlinear  range.  These  instabilities  are  characterized  by  an  unbounded  growth  of 
the  energy  in  finite  time,  despite  the  dissipative  character  of  the  underlying  physical  sys¬ 
tems  under  consideration.  This  situation  clearly  identifies  the  need  for  new  and  improved 
numerical  algorithms,  motivating  the  developments  considered  in  this  project. 

The  completely  different  nature  of  inelastic  problems,  involving  the  additional  set  of  plas¬ 
tic/damage  evolution  equations  (usually  of  a  unilaterally  constrained  character  due  to  the 
presence  of  the  so-called  yield/damage  surface  condition),  required  the  development  of  new 
ideas  and  schemes  beyond  the  ones  available  for  the  elastic  case.  Our  approach  in  address¬ 
ing  these  challenges  can  be  summarized  in  the  actual  analysis  of  the  discrete  equations 
driving  the  design  of  the  new  numerical  schemes,  and  the  exact  incorporation  (by  design) 
in  the  discrete  system  of  the  energy- moment  um  response  of  the  underlying  physical  system. 
In  this  way,  the  rigorous  proof  of  all  the  conservation/dissipation  properties  for  general 
material  models,  accounting  also  for  the  finite  element  interpolation  of  the  governing  equa¬ 
tions,  provides  the  robustness  required  to  the  computational  tools  for  the  analysis  of  the 
complex  practical  applications  of  interest  to  the  Air  Force.  In  this  respect,  we  note  that 
the  mathematical  analyses  of  the  algorithms  identified  the  need  to  develop  new  assumed 
strain  finite  element  methods  that  conformed  with  the  dissipative/conservation  properties 
of  the  temporal  approximation.  This  situation  identified  a  new  objective  of  the  project, 
also  successfully  achieved. 

The  improved  stability  properties  of  the  newly  developed  integration  algorithms  are  then 
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supported  by  rigorous  mathematical  analyses  of  the  discrete  dynamical  systems  that  they 
generate.  We  also  have  given  a  special  attention  to  the  actual  implementation  of  the  new 
algorithms  in  the  context  of  the  finite  element  method.  The  combination  of  all  these  results 
has  led  to  powerful  novel  computational  tools,  with  the  sound  theoretical  basis  necessary 
for  the  analysis  of  the  complex  practical  problems  of  interest  to  the  Air  Force. 


2.  Brief  Description  of  the  Accomplishments 

The  major  results  accomplished  under  this  project  can  be  summarized  as  follows: 

1.  Energy-dissipative,  momentum- conserving  algorithms  (EDMC  schemes)  for  finite  strain 
multiplicative  plasticity.  Continuing  with  preliminary  results  of  a  previous  AFOSR  research 
effort,  we  have  developed  new  integration  algorithms  for  finite  strain  plasticity  that  incor¬ 
porate  the  characteristic  (and  critical)  property  of  non-negative  energy  dissipation  in  the 
discrete  system,  while  inheriting  also  the  conservation  laws  of  linear  and  angular  momenta. 
The  construction  of  the  new  algorithms  makes  a  crucial  use  of  the  arguments  employed 
in  the  characterization  of  the  the  physical  and  mathematical  properties  of  the  continuum 
model,  in  the  context  of  modern  treatments  based  on  a  multiplicative  decomposition  of 
the  deformation  gradient  in  an  elastic  and  plastic  part.  The  resulting  scheme  conserves 
exactly  both  the  linear  and  angular  momenta,  and  provides  the  exact  physical  dissipation 
predicted  by  the  model,  including  the  exact  energy  conservation  for  an  elastic  step.  This 
leads  not  only  to  physically  better  numerical  results,  but  also  to  more  numerically  stable 
simulations.  This  improved  numerical  stability  is  to  be  contrasted  with  the  numerical  in¬ 
stabilities  observed  by  existing  schemes  in  this  strongly  nonlinear  problem,  even  with  its 
physically  dissipative  character.  These  results  appeared  in  [1,2]  1 . 

2.  Volume-preserving  energy-momentum  schemes  for  isochoric  plastic  models.  Several  im¬ 
portant  applications  of  elastoplastic  models  (namely,  metals)  exhibit  an  isochoric  plastic 
response,  that  is,  there  is  no  plastic  change  of  volume.  This  response  is  a  very  characteristic 
physical  property  of  these  materials  and,  as  such,  crucial  to  reproduce  in  the  actual  numer¬ 
ical  simulations.  This  physical  property  translates  in  a  very  rich  geometric  structure  of  the 
aforementioned  multiplicative  elastoplastic  models,  with  a  key  role  played  by  the  elastic 
metric  in  the  intermediate  configuration  defined  by  the  plastic  part  of  the  deformation  gra¬ 
dient.  We  have  developed  a  new  time-stepping  approximation  of  this  geometric  structure 
that,  when  combined  with  the  developments  summarized  above  for  general  elastoplastic 
models,  results  in  a  volume  preserving  scheme  that  also  exhibits  the  exact  non-negative 
energy  dissipation  and  momentum  conservation  laws  of  the  underlying  physical  system. 
The  improved  numerical  stability  comes  then  with  the  added  benefit  of  reproducing  the 
plastic  isochoric  response  of  interest.  The  new  algorithm  defines  a  new  approximation 

t  The  numbering  of  the  references  follows  the  list  of  publications  presented  below  in 
page  14  of  this  report. 
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FIGURE  2.1  Three-dimensional  elastoplastic  solid  in  free  flight.  Sequence  of  deformations 
in  the  early  stages  of  the  motion  with  the  spatial  distribution  of  the  equivalent  plastic  strain. 
Solution  obtained  with  the  new  EDMC  scheme. 
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FIGURE  2.2  Three-dimensional  elastoplastic  solid  in  free  flight.  Temporal  evolutions  of  the 
the  three  components  of  the  linear  momentum  (left)  and  the  angular  momentum  (right)  for  the 
solution  obtained  with  the  new  EDMC  scheme,  showing  the  exact  conservation  of  all  these  quan¬ 
tities.  They  correspond  exactly  to  the  momenta  associated  with  the  assumed  initial  velocities 
for  the  assumed  undeformed  initial  configuration. 


TRAPEZOIDAL  RULE 


EDMC  SCHEME 


Time 


Time 


FIGURE  2.3  Three-dimensional  elastoplastic  solid  in  free  flight.  Temporal  evolutions  of  the 
total  energy  obtained  by  the  trapezoidal  rule  (left)  and  EDMC  scheme  (right)  for  different 
time  step.  The  non-negative  energy  dissipative  character  of  the  new  EDMC  scheme  is  to  be 
contrasted  with  the  persistent  energy  growth  observed  with  the  trapezoidal  rule  leading  to  the 
termination  of  the  simulation  for  a  lack  of  convergence  of  the  Newton-Raphson  scheme  used  in 
the  solution  process. 
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of  the  aforementioned  geometric  structure  that  controls  the  volumetric  strain.  The  new 
algorithm  is  second  order  accurate  in  time,  as  desired. 

Figures  2.1,  2.2  and  2.3  illustrate  these  developments.  They  show  the  solution  obtained 
in  the  modeling  of  the  free  flight  of  a  three-dimensional  solid  consisting  of  a  very  stiff 
elastic  cylindrical  core  and  two  flexible  arms  made  of  an  elastoplastic  material  modeled 
with  J2-flow  theory  of  finite  strain  plasticity  following  the  classical  Mises  yield  function. 
The  elastic  response  is  based  on  Hencky’s  hyperelastic  potential  on  the  logarithmic  elastic 
strains.  The  change  of  volume  is  then  entirely  elastic,  with  no  plastic  volumetric  strain. 
The  solid  is  given  an  initial  velocity  corresponding  to  a  rotation  around  its  center.  The 
solid  is  in  free  flight  afterwards.  The  total  linear  and  angular  momenta  are  then  exactly 
conserved  in  this  free  motion,  with  the  energy  being  dissipated  (reduced)  if  plastic  flow 
occurs  while  being  fully  conserved  during  purely  elastic  evolution  steps.  Figure  2.1  shows 
the  deformed  configuration  of  the  solid  at  different  times  computed  with  the  new  volume¬ 
preserving  EDMC  scheme,  depicting  also  the  distribution  of  the  equivalent  plastic  strain 
(thus  showing  the  extent  of  the  plastic  flow  in  the  material).  The  large  finite  deformations 
are  apparent.  Figure  2.2  shows  the  evolution  in  time  of  the  three  components  of  the 
linear  momentum  (left)  and  angular  momentum  (right)  of  the  solid,  confirming  their  exact 
conservation.  Figure  2.3  shows  the  evolution  of  the  total  energy  of  the  solution  computed 
with  the  EDMC  scheme  and  the  classical  trapezoidal  rule  (a  member  of  the  classical 
family  of  Newmark  schemes)  in  combination  with  an  exponential  approximation  of  the 
plastic  evolution  equations,  all  for  different  time  steps.  A  non-increasing  energy  evolution 
should  be  obtained  after  the  initial  stages  when  the  load  is  applied.  We  can  observe 
that,  despite  the  dissipative  character  of  the  underlying  physical  system,  the  trapezoidal 
rule  shows  an  uncontrolled  growth  of  energy  (blow-up)  forcing  to  stop  the  computation. 
This  numerical  instability  is  observed  for  different  time  step  sizes.  This  situation  is  to  be 
contrasted  with  the  energy  evolution  for  the  solution  obtained  by  the  new  EDMC  scheme. 
The  energy  dissipation  is  always  strictly  non- negative  (by  design)  and,  in  fact,  exact.  In 
particular,  full  energy  conservation  is  obtained  during  the  elastic  steps  (again,  by  design). 
The  significantly  improved  numerical  stability  properties  of  the  new  scheme  are  apparent. 
Furthermore,  plots  (not  shown)  confirm  the  lack  of  any  plastic  change  of  volume  in  the 
deformation  of  the  solid,  thus  exactly  reproducing  the  response  of  the  solid. 

We  presented  the  development  of  the  volume-preserving  energy-dissipative  momentum- 
conserving  algorithms  for  isochoric  models  of  finite  strain  plasticity  in  the  papers  [3,5].  A 
summary  of  these  results,  with  additional  details  of  the  ones  discussed  above,  can  be  found 
in  Appendix  I  of  this  final  report. 

3.  Conserving  assumed  strain  finite  elements  for  continuum,  problems.  The  quasi-incompre- 
ssible  character  of  the  plastic  response  of  metals  requires  the  consideration  of  locking-free 
finite  elements  for  their  spatial  discretization.  Assumed  strain  methods  treating  the  volu¬ 
metric  strain  in  a  separate  manner  have  been  proven  to  be  a  reliable  option  in  these  situa¬ 
tions.  However,  the  direct  use  of  existing  techniques  destroys  the  conservation/dissipation 
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FIGURE  2.4  Impact  of  a  cylindrical  copper  rod  on  a  rigid  wall,  a)  Problem  definition  (dis¬ 
tances  in  mm),  showing  a  quarter  of  the  specimen  discretized  with  976  brick  elements.  The  rod 
is  given  the  initial  axial  velocity  of  v0  =  0.227  mm/ /is ,  modeling  the  impact  on  the  rigid  wall 
by  imposed  axial  displacements  at  its  base.  b)  Evolution  of  the  energy  (total  with  potential 
and  kinetic)  showing  the  monotonic  decrease  of  the  total  energy  consequence  of  the  exactly  cap¬ 
tured  plastic  dissipation  by  the  new  EDMC  time-stepping  scheme  in  combination  with  the  new 
assumed  strain  finite  elements.  Existing  methods  require  a  large  amount  of  artificial  numerical 
dissipation  due  to  the  appearance  of  highly  oscillatory  and  unstable  solutions,  c)  Deformed 
configurations  in  time  showing  the  distribution  of  the  equivalent  plastic  strain  (half  the  speci¬ 
men  shown  from  below- front).  The  large  strains  are  apparent  (illustrating  the  lack  of  locking), 
and  so  is  the  bulging  of  the  specimen  as  the  specimen  shortens,  agreeing  with  experimental 
tests  for  the  considered  material  of  pure  copper  with  large  strain  hardening. 
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properties  of  the  temporal  discretizations  discussed  above,  thus  requiring  the  development 
of  new  alternative  techniques.  We  have  developed  new  assumed  strain  finite  elements 
based  on  the  classical  scaling  of  the  volumetric  part  of  the  deformation  gradient  (i.e.  its 
determinant  or  Jacobian)  by  a  mixed  Jacobian,  but  requiring  an  alternative  definition 
of  the  associated  linearized  strain  operator  in  the  equations  of  motion.  After  identify¬ 
ing  the  conditions  for  the  final  scheme  to  be  numerically  consistent  and  to  comply  with 
the  conservation  laws  of  linear  and  angular  momenta,  the  associated  relative  equilibria, 
and  the  evolution  of  the  energy,  a  general  strategy  for  the  definition  of  this  “conserving 
consistent  linearization”  has  been  identified.  Based  on  this  new  approach,  we  have  devel¬ 
oped  new  quadrilateral  and  triangular  elements  for  plane  problems,  and  brick  elements  for 
three-dimensional  applications.  For  example,  a  fully  EDMC  trilinear  brick  with  constant 
volumetric  strain  (Q1/A0)  and  quadratic  quad  with  linear  volumetric  strain  (Q2/A1)  are 
now  available. 

Figure  2.4  illustrates  the  performance  of  the  newly  developed  conserving  assumed  strain 
finite  elements  in  a  typical  application  with  large  (isochoric)  plastic  strains.  It  shows  the 
results  obtained  in  the  classical  benchmark  problem  of  the  impact  of  a  cylindrical  rod 
on  a  rigid  wall  (Taylor’s  problem).  The  new  EDMC  time-stepping  algorithms  avoid  any 
instabilities  in  time,  in  contrast  with  classical  existing  schemes  that  require  an  excessive 
amount  of  numerical  dissipation  to  accomplish  a  meaningful  solution  (without  spurious 
oscillations  and  instabilities).  We  note  again  that  this  is  a  direct  consequence  of  the  new 
schemes  capturing  exactly  the  physical  dissipation  occurring  in  the  solid.  The  effective¬ 
ness  of  the  new  assumed  strain  finite  elements  in  avoiding  volumetric  locking  is  confirmed 
by  the  large  strains  involved  and  their  agreement  with  experimental  observations  (not 
shown).  Basic  displacement  finite  element  methods  would  lock  in  this  simulation  due  to 
the  quasi-incompressible  character  of  the  deformation  given  the  isochoric  plastic  response 
of  the  solid.  The  new  elements  are  locking  free  while  exactly  exhibiting  the  energy  dissi¬ 
pation  and  momentum  conservation  properties  of  the  EDMC  scheme  used  in  the  temporal 
discretization. 

These  results  were  presented  in  the  papers  [7,8,10].  Additional  details  can  be  found  in 
Appendix  II  of  this  final  report. 

4.  Energy-dissipative  momentum-conserving  algorithms  for  finite  strain  thermo- elasticity. 
We  have  developed  new  time-stepping  algorithms  for  coupled  thermo-elasticity  that  in¬ 
herit  by  design  the  a-priori  stability  estimates  of  this  physical  system  and  its  momenta 
conservation  laws.  The  stress  formula  developed  previously  for  uncoupled  cases  has  been 
extended  to  account  for  the  temperature  coupling,  defining  in  the  process  the  proper  tem¬ 
poral  approximation  of  the  conjugate  entropy  field.  The  proper  treatment  of  the  energy 
evolution  equation,  including  heat  conduction,  has  allowed  to  arrive  to  a  new  numerical 
scheme  that  exhibits  rigorously  the  dissipative  character  of  the  so-called  canonical  free  en¬ 
ergy  characteristic  of  these  systems.  Complete  mathematical  analyses  are  available.  These 
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FIGURE  2.5  Tumbling  of  a  L-shaped  thermo-elastic  block:  solution  obtained  with  the  new 
EDMC  schemes  for  thermo-elasticity,  a)  Motion  depicted  by  the  configuration  of  the  solid  at 
different  times  (left  to  right  and  down)  with  the  distribution  of  the  temperature.  The  unde¬ 
formed  solid  (first  top  left  frame)  is  given  an  initial  loading  on  the  top  and  bottom  faces,  being 
in  free  flight  after  this  initial  phase.  The  series  of  deformations  illustrate  the  large  displacements 
and  strains  that  the  solid  is  subjected  to,  with  the  temperature  variation  occurring  from  the 
associated  thermo-elastic  coupling,  b)  Evolution  of  the  three  components  of  the  angular  mo¬ 
mentum.  After  the  initial  loading  phase,  the  three  components  are  exactly  conserved,  as  they 
must  physically,  showing  the  momentum  conserving  properties  of  the  new  scheme,  c)  Evolution 
of  the  canonical  free  energy  with  time.  The  new  EDMC  scheme  conforms  with  the  a-priori  sta¬ 
bility  estimate  defined  by  the  monotonic  decrease  of  the  canonical  free  energy,  in  contrast  to  the 
solution  obtained  by  a  classical  scheme  like  trapezoidal  rule  exhibiting  a  nonlinear  numerical 
instabilities  in  the  form  of  an  unbounded  growth  of  the  energy. 
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schemes  can  also  be  incorporated  in  a  staggered  (partitioned)  solution  of  the  problem. 
These  results  have  been  compiled  in  [12],  with  the  plans  of  additional  publications.  In 
particular,  the  extension  to  coupled  thermo-plasticity  is  also  under  consideration. 

Figure  2.5  illustrates  the  performance  of  the  new  EDMC  scheme  developed  in  this  part 
of  the  project.  It  includes  the  deformation  and  temperature  distribution  of  a  L-shaped 
thermo-elastic  block  tumbling  in  space.  The  conservation  of  the  angular  momentum  and 
the  monotonic  decrease  of  the  canonical  free  energy  due  to  heat  conduction  after  the  initial 
loading  conforms  with  the  same  properties  of  the  continuum  system.  Thermal  heat  con¬ 
duction  introduces  dissipation  in  the  physical  system,  non-negative  from  the  fundamental 
second  law  of  thermodynamics,  and  leading  to  the  monotonic  decrease  of  the  so-called 
canonical  free  energy.  This  defines  an  a-priori  stability  estimate  that  the  numerical  simu¬ 
lation  must  conform  to  for  numerical  stability  (in  the  linear  case,  the  canonical  free  energy 
strictly  defines  a  norm  of  the  solution  in  the  displacement,  velocity  and  temperature).  The 
plot  shows  that  the  new  EDMC  does  conform  with  this  fundamental  energy  principle, 
thus  leading  to  stable  simulations  of  the  motions  of  interest.  This  property  is  supported 
by  rigorous  mathematical  analyses  rigorously  showing  its  satisfaction.  This  situation  is 
to  be  contrasted  with  the  performance  of  standard  numerical  schemes  like  the  classical 
trapezoidal  rule,  also  shown  in  the  figure.  The  energy  evolution  is  not  monotonic  in  this 
case  and,  in  fact,  it  leads  to  an  unbounded  growth,  thus  showing  the  instability  of  these 
classical  schemes  in  the  considered  nonlinear  range.  The  simulation  blows  up  in  this  case, 
a  performance  that  needs  to  be  contrasted  again  with  the  new  EDMC  scheme.  These 
instabilities  are  completely  and  rigorously  avoided  in  the  new  schemes  developed  in  this 
project. 

5.  Extension  to  shell  models.  We  have  also  considered  the  application  of  the  newly  de¬ 
veloped  time-stepping  algorithms  to  problems  involving  shells,  as  this  is  one  of  the  main 
configurations  in  Air  Force  applications.  In  this  way,  we  have  developed  the  implemen¬ 
tation  of  these  new  methods  in  reduced-solid  shell  elements,  incorporating  also  the  new 
assumed  strain  treatment  developed  for  continuum  problems  discussed  above  (Item  3). 
Enhanced  strain  treatments  for  the  through-the-thickness  shell  response  as  well  as  the 
incorporation  of  common  assumed  strain  treatments  of  the  bending/shear  response  have 
been  considered.  We  have  considered  both  plastic  and  thermoelastic  shells.  In  both  cases, 
the  implementation  follows  the  ideas  presented  in  [8]  for  the  construction  of  assumed  strain 
finite  element  methods. 

Figure  2.6  illustrate  these  results,  currently  under  further  development.  It  shows  the 
solution  obtained  in  the  modeling  of  the  motion  of  a  L-shaped  elastoplastic  shell.  The 
shell  is  loaded  by  transversal  tractions  along  its  extreme  edges  and  at  the  central  corner 
for  a  period  of  time,  being  in  free-flight  after  this  initial  loading  phase.  The  result  ing 
motion  is  shown  in  that  figure  depicting  the  deformed  configuration  of  the  shell  at  different 
times.  The  large  displacement,  rotations  and  strains  are  apparent.  The  distribution  of 


Final  Report.  FA9550-05-1-0117 


11 


a) 


b) 

Angular  Momentum 


c) 

Energy 


FIGURE  2.6  Tumbling  of  a  L-shaped  elastoplastic  shell:  solution  obtained  with  the  new 
EDMC  schemes  for  finite  strain  plasticity  of  nonlinear  shells,  a)  Motion  depicted  by  the  con¬ 
figuration  of  the  solid  at  different  times  (left  to  right  and  down)  with  the  distribution  of  the 
equivalent  plastic  strain.  The  undeformed  solid  (first  top  left  frame)  is  loaded  with  distributed 
forces  at  both  extreme  edges  and  central  corner,  being  in  free-flight  after  the  initial  loading 
phase.  The  series  of  deformations  illustrate  the  large  displacements  and  strains  that  the  solid 
is  subjected  to,  including  large  plastic  strains,  b)  Evolution  of  the  three  components  of  the 
angular  momentum.  After  the  initial  loading  phase,  the  three  components  are  exactly  conserved 
as  it  must  happen  from  physical  considerations,  hence  showing  the  momentum  conserving  prop¬ 
erties  of  the  new  scheme,  c)  Evolution  of  the  total  energy  and  its  kinetic  and  potential  (strain 
plus  hardening)  energies.  The  monotonic  evolution  of  the  total  energy  is  to  be  noted,  with  the 
energy  decrease  corresponding  to  the  plastic  dissipation  captured  exactly  by  the  new  EDMC 
scheme. 
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the  equivalent  plastic  strain  is  also  included,  showing  the  extent  of  the  plasticity.  The 
important  aspect  is  that  the  new  EDMC  scheme  is  able  to  capture  exactly  the  energy 
dissipation  of  this  plastic  response.  The  evolution  of  the  energy  included  in  Figure  2.6 
confirms  this  property. 

6.  Other  results:  numerical  modeling  of  fracture  and  failure  through  the  strong  discon¬ 
tinuity  approach.  The  previous  accomplishments  defined  the  main  focus  of  this  research 
project  as  outlined  in  the  original  grant  proposal.  In  addition,  we  have  started  exploring 
the  use  of  finite  elements  with  embedded  strong  discontinuities  for  the  analysis  of  dynamic 
failure  of  materials.  The  approach  enhances  the  interpolations  of  traditional  finite  elements 
with  the  discontinuous  solutions  characteristic  of  the  ultimate  stages  of  the  deformation  of 
solids  (e.g.  cracks  or  shear  bands).  A  multiscale  approach  allows  to  introduce  these  solu¬ 
tions  (also  known  as  strong  discontinuities)  into  the  finite  elements  at  the  local  level,  with 
enhanced  parameters  eliminated  locally  through  their  static  condensation  hence  resulting 
in  a  very  efficient  solution  of  the  global  mechanical  problem.  A  main  advantage  is  that  no 
remeshing  is  needed  for  the  resolution  of  the  discontinuity. 

We  have  presented  in  [4,6]  a  new  strategy  for  the  development  of  these  enhancements  in 
the  infinitesimal  range  of  small  strains.  The  approach  is  based  on  the  direct  introduction  in 
the  discrete  strain  field  of  the  finite  elements  of  the  strain  modes  necessary  for  the  accurate 
resolution  of  the  kinematics  associated  with  a  discontinuous  displacement  held.  Piece-wise 
linear  interpolations  of  the  discontinuity  jumps  are  considered  along  the  discontinuity,  in 
contrast  with  the  piece-wise  constant  interpolations  existing  to  date.  The  new  numer¬ 
ical  approach  leads,  in  particular,  to  quadrilateral  finite  elements  able  to  represent  the 
separation  of  the  discontinuities  without  a  spurious  transfer  of  stresses  (or  stress  locking). 

The  extension  of  these  considerations  to  the  finite  deformation  has  been  presented  in 
[9].  In  this  case,  the  key  aspect  has  been  the  development  of  a  new  enhancement  of  the 
deformation  gradient,  that  conforms  with  the  geometric  structure  of  the  problem,  leading 
in  particular  to  frame  indifferent  formulations.  Furthermore,  the  new  strategy  has  shown  to 
be  particularly  appropriate  for  the  analysis  of  failure  and  fracture  of  solids  in  the  dynamic 
range.  Preliminary  results  can  be  found  in  the  technical  report  [11]  for  the  infinitesimal 
case.  These  results  have  motivated  and  defined  the  research  lines  to  follow  in  this  effort  as 
describe  in  the  following  section. 


3.  Current  and  Future  Work 

The  results  obtained  in  this  project  have  allowed  to  identify  a  number  of  extensions  and 
new  lines  of  research  for  the  near  future.  As  indicated  at  the  end  of  the  last  section,  the 
numerical  analysis  of  the  failure  and  fracture  of  solids  defines  an  important  problem  to 
be  addressed,  given  the  many  applications  of  interest  to  the  Air  Force.  In  particular,  the 
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numerical  resolution  of  the  highly  non-smooth  solutions  observed  in  the  ultimate  stages 
of  the  deformation  of  solids  sets  out  a  very  challenging  problem.  These  solutions  involve 
strong  discontinuities  whose  modeling  in  a  numerically  discrete  context  require  the  design 
of  special  numerical  methods  for  their  resolution.  Typical  examples  of  these  considerations 
are  cracks  in  brittle  or  quasi-brittle  materials,  or  shear  bands  (and  similar  localization 
bands)  in  ductile  failures. 

As  noted  above,  we  are  already  involved  in  the  development  of  new  finite  elements  to 
capture  these  discontinuities  at  the  element  interiors.  This  is  accomplished  through  a 
multi-scale  treatment  of  the  mechanical  problem  at  hand,  with  the  consideration  of  the 
discontinuities  at  the  local  (element)  level,  thus  leading  to  very  efficient  numerical  tech¬ 
niques,  easily  accommodating  existing  numerical  methods  for  the  large  scale. 

The  dynamic  range  defines  a  particularly  difficult  setting  for  the  different  additional  effects 
to  be  resolved  accurately  and  in  a  stable  manner.  For  example,  crack  branching  defines 
a  difficult  problem  to  be  captured  numerically,  because  of  both  the  spatial  discretization 
and  the  time-stepping  scheme.  In  the  former,  and  in  the  context  of  the  considered  finite 
elements  with  embedded  strong  discontinuities,  the  design  of  new  element  enhancements 
that  can  accommodate  a  bifurcating  discontinuity  defines  a  clear  objective  for  the  con¬ 
tinuing  effort.  Similarly,  the  development  of  stable  temporal  discretizations  (extending 
the  energy-dissipative  momentum-conserving  schemes  developed  in  this  project  to  cases 
involving  an  inelastic  cohesive  law  along  discontinuity  surface)  is  another  clear  objective 
for  future  work.  In  addition,  the  need  to  consider  additional  physical  effects  (like  coupled 
thermoplastic  responses  at  high  strain  rates)  motivates  also  the  consideration  of  the  ex¬ 
tension  of  the  work  develop  here  to  those  cases  (e.g.  EDMC  time-stepping  algorithms  for 
nonlinear  coupled  thermoplasticity). 

This  continuing  effort  will  widen  the  range  of  application  to  the  methods  that  we  have 
developed  in  the  current  research  project.  The  overall  approach  followed  in  it  (that  is, 
the  mathematical  analysis  of  the  discrete  equations  drives  the  design  of  new  numerical 
schemes)  will  be  crucial  in  all  of  our  future  developments  in  this  area. 
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seminar,  Department  of  Aerospace  and  Mechanical  Engineering,  University  of  South¬ 
ern  California,  January  20,  2005. 

2.  “Energy-Dissipative  Momentum-Conserving  Time-Stepping  Algorithms  for  Dynamic 
Plastic  Finite  Strain  Plasticity,”  invited  contribution,  VIII  International  Conference 
on  Computational  Plasticity  (COMPLAS  8),  Barcelona,  Spain,  September  4-6  2005. 

3.  “Finite  Element  Locking  and  the  Enhanced  Strain  Formulation,”  CISM  course  on 
Mixed  Finite  Element  Technologies,  CISM  Udine,  Italy,  September  26-30,  2005. 

4.  “Numerical  Integration  of  the  Nonlinear  Dynamics  of  Inelastic  Solids  and  Structures,” 
Department  of  Civil  and  Environmental  Engineering,  University  of  California  at  Los 
Angeles  (UCLA),  February  14  2006. 

5.  “Numerical  Integration  of  the  Nonlinear  Dynamics  of  Elastoplastic  Solids,”  keynote 
lecture,  3rd  European  Conference  on  Computational  Mechanics  (ECCM  3),  Lisbon, 
Portugal,  June  5-9  2006. 

6.  “Energy-Momentum  Schemes  for  Finite  Strain  Plasticity,”  keynote  lecture,  7th  World 
Congress  on  Computational  Mechanics  (WCCM  7),  Los  Angeles  CA,  July  17-21  2006. 

7.  “Finite  Elements  with  Embedded  Strong  Discontinuities  of  Higher  Order  Kinematics,” 
invited  contribution,  7th  World  Congress  on  Computational  Mechanics  (WCCM  7), 
Los  Angeles  CA,  July  17-21  2006. 

8.  “Recent  Developments  in  the  Formulation  of  Finite  Elements  with  Embedded  Strong 
Discontinuities,”  invited  contribution,  IUTAM  Symposium  on  Discretization  Methods 
for  Evolving  Discontinuities,  Lyon,  France,  September  4-7  2006. 

9.  “Energy-Momentum  Algorithms  for  Nonlinear  Solid  Dynamics  and  their  Assumed 
Strain  Finite  Element  Formulation,”  invited  semi-plenary  lecture,  Computational  Metli- 
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ods  in  Structural  Dynamics  and  Earthquake  Engineering  (COMPDYN07),  Rethymno, 
Greece,  June  13-16  2007. 

10.  “New  Finite  Elements  with  Embedded  Strong  Discontinuities,”  invited  contribution, 
9th  US  National  Congress  on  Computational  Mechanics  (USNCCM  9),  San  Francisco, 
CA,  July  23-26,  2007. 

11.  “Dynamic  Fracture  Using  Finite  Elements  Enhanced  with  Cohesive  Discontinuities,” 
9th  US  National  Congress  on  Computational  Mechanics  (USNCCM  9),  San  Francisco, 
CA,  July  23-26,  2007. 

12.  “Finite  Elements  with  Embedded  Discontinuities  and  Dynamic  Fracture,”  invited  con¬ 
tribution,  IX  International  Conference  on  Computational  Plasticity  (COMPLAS  9), 
Barcelona,  Spain,  September  5-7  2007. 

13.  “Numerical  Integration  in  Nonlinear  Solid  and  Structural  Dynamics,”  Department  of 
Civil  and  Environmental  Engineering,  Johns  Hopkins  University,  December  17  2007. 

14.  “Energy-Momentum  Algorithms  for  Nonlinear  Dynamics  of  Elastoplastic  Solids,”  in¬ 
vited  contribution,  IUTAM  Symposium  on  Theoretical,  Modeling  and  Computational 
Aspects  of  Inelastic  Media,  Cape  Town,  South  Africa,  January  14-18  2008. 

15.  “Finite  Element  Modeling  of  Kirchlioff  Rods,”  invited  contribution,  6th  International 
Conference  on  Computation  of  Shell  and  Spatial  Structures  (IASS-IACM  2008),  Cor¬ 
nell  University,  Ithaca,  NY,  May  28-31  2008. 

16.  “Energy-Momentum  Algorithms  for  Nonlinear  Coupled  Thermo-Elastodynamics,”  in¬ 
vited  keynote  lecture,  8th  World  Congress  on  Computational  Mechanics  (WCCM  8), 
Venice,  Italy,  June  30-July  4  2008. 

17.  “Modeling  of  Dynamic  Fracture  using  Finite  Elements  with  Embedded  Strong  Dis¬ 
continuities,”  invited  contribution,  8th  World  Congress  on  Computational  Mechanics 
(WCCM  8),  Venice,  Italy,  June  30-July  4  2008. 


7.  Honors  and  Awards 

The  PI,  Francisco  Armero  was  elected  Fellow  of  the  International  Association  of  Compu¬ 
tational  Mechanics  (IACM)  in  2004.  The  Fellows  Award  was  given  to  him  during  the  6th 
World  Congress  on  Computational  Mechanics  that  took  place  in  Beijing,  China,  September 
5-11,  2004.  Additionally,  Francisco  Armero  was  awarded  in  the  past  the  Young  Investi¬ 
gator  Award  by  the  International  Association  for  Computational  Mechanics  (IACM)  in 
July  2002,  the  Juan  C.  Simo  award  and  medal  by  SEMNI  (Spanish  Society  for  Numer¬ 
ical  Methods  in  Engineering)  in  June  1999,  the  NSF  CAREER  award  (National  Science 
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Foundation)  in  June  1997,  and  the  ONR  Young  Investigator  Award  (Office  of  the  Naval 
Research)  in  June  1996.  He  also  received  the  best  paper  award  for  “the  most  outstanding 
paper  published  in  Engineering  Computations  in  the  year  1997” . 

Additional  honors  to  the  PI,  Francisco  Armero,  related  to  this  research  project  during  the 
performance  of  the  grant  include  service  as: 

1.  Member  of  the  editorial/advisory  board  of: 

•  Communications  in  Numerical  Methods  in  Engineering,  April  2005-present. 

•  ASCE  Journal  of  Engineering  Mechanics,  Associate  Editor,  Sept.  2003-Oct.  2005. 

•  Informes  de  la  Construccion,  December  2006-present. 

•  Computers  &  Concrete,  April  2003-present. 

•  Finite  Elements  in  Analysis  and  Design,  February  2002-present. 

•  Computer  Methods  in  Applied  Mechanics  and  Engineering,  June  2001-present. 

•  International  Journal  for  Numerical  Methods  in  Engineering,  January  2001-present. 

•  Computers  &  Structures,  November  1998-present. 

•  International  Journal  of  Numerical  Methods  in  Fluids,  November  1997-January 
2008. 

2.  Committee  service: 

•  Committee  on  Computational  Mechanics,  ASCE  Engineering  Mechanics  Division 
(August  1999-present),  vice-chair  (2003-2004,  2006- present),  chair  (August  2004- 
2006). 

•  Technical  Advisory  Committee,  10th  International  Conference  on  Computational 
Plasticity  (COMPLAS  X),  Barcelona,  Spain,  September  2-4,  2009. 

•  International  Advisory  Board,  Computational  Methods  in  Structural  Dynamics  and 
Earthquake  Engineering  (COMPDYN09),  Island  of  Rhodes,  Greece,  June  22-24, 
2009. 

•  International  Advisory  Board,  1st  International  Conference  on  Computational  Tech¬ 
nologies  in  Concrete  Structures  (CTCS09),  Seoul,  Korea,  May  2009. 

•  Scientific  Advisory  Committee,  9th  International  Conference  on  Computational 
Structures  Technology,  Athens,  Greece,  September  2-5,  2008. 

•  Local  Organizing  Committee,  9  US  National  Congress  on  Computational  Mechanics 
(USNCCM  9),  San  Francisco  CA,  July  17-19,  2007. 

•  International  Advisory  Board,  Computational  Methods  in  Structural  Dynamics  and 
Earthquake  Engineering  (COMPDYN07),  Rethymnon,  Crete,  Greece,  June  13-15, 
2007. 

•  Technical  Advisory  Committee,  International  Conference  on  Coupled  Problems, 
Ibiza,  Spain,  May  21-23,  2007. 
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•  Technical  Advisory  Committee,  9th  International  Conference  on  Computational 
Plasticity  (COMPLAS  IX),  Barcelona,  Spain,  September  5-7,  2007. 

•  Scientific  Programme  Committee,  7th  World  Congress  on  Computational  Mechanics 
(WCCM7),  Los  Angeles  CA,  U.S.A.,  July  16-22,  2006. 

•  Scientific  Committee,  III  European  Conference  on  Computational  Mechanics,  Lis¬ 
bon,  Portugal,  June  5-9,  2006. 

•  Scientific  Advisory  Committee,  8th  International  Conference  on  Computational  Struc¬ 
tures  Technology,  Las  Palmas  de  Gran  Canaria,  Spain,  September  12-15,  2006. 

•  Scientific  Committee,  5th  International  Conference  on  Computation  of  Shells  and 
Spatial  Structures  (IASS-IACM  2005),  Salzburg,  Austria,  June  1-4,  2005. 

•  Technical  Advisory  Committee,  International  Conference  on  Coupled  Problems, 
Santorini  Island,  Greece,  May  25-28,  2005. 

•  Technical  Advisory  Committee,  8th  International  Conference  on  Computational 
Plasticity  (COMPLAS  VIII),  Barcelona,  Spain,  September  5-8,  2005. 

•  Scientific  Advisory  Committee,  7th  International  Conference  on  Computational  Struc¬ 
tures  Technology,  Lisbon,  Portugal,  September  7-9,  2004. 

3.  Organizer  of  symposia: 

•  “Numerical  Techniques  for  the  Modeling  of  Material  Failure,”  7  sessions,  38  contri¬ 
butions,  8th  World  Congress  on  Computational  Mechanics  (VIII  WCCM),  Venice, 
Italy,  June  30-July  4,  2008. 

•  “Numerical  Techniques  for  the  Modeling  of  Material  Failure  in  Solids:  Symposium 
in  Honor  of  Professor  Kaspar  Wiliam  on  the  Occasion  of  his  65th  birthday,”  7  ses¬ 
sions,  32  contributions,  9th  US  National  Congress  on  Computational  Mechanics  (IX 
USNCCM),  San  Francisco  CA,  July  23-26  2007. 

•  “Modeling  and  Numerical  Simulation  of  of  Failure  in  Inelastic  Shells  and  Spatial 
Structures,”  5th  International  Conference  on  Computation  of  Shell  and  Spatial 
Structures  (IASS-IACM  2005),  Salzburg,  Austria,  June  1-4  2005. 


8.  Outline  of  the  Rest  the  Report 

We  have  added  two  appendices  to  this  report  to  illustrate  in  more  detail  some  of  the 
technical  results  obtained  in  this  project.  Appendix  I  describes  the  formulation  of  energy- 
dissipative  momentum-conserving  schemes  for  finite  strain  plasticity  that  also  conserve 
the  plastic  volumetric  response  of  the  underlying  model  (Item  2  in  Section  2).  Appendix 
II  describes  the  development  of  the  new  assumed  strain  finite  element  methods  for  the 
conserving  time-stepping  algorithms  (Item  3  in  Section  2).  A  more  detailed  description 
follows. 
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8.1.  Appendix  I:  Volume-Preserving  Energy-Dissipative  Momentum- 
Conserving  Time-Stepping  Algorithms  for  Isochoric  Plastic  Models 
of  Multiplicative  Plasticity 

This  appendix  presents  a  new  energy-dissipative  momentum- conserving  algorithm  for  mul¬ 
tiplicative  finite  strain  plasticity  (F  =  FeFp)  that  also  preserves  exactly  the  plastic  volume 
for  isochoric  plastic  models.  The  new  algorithm  exhibits  exactly  the  conservation  laws  of 
linear  and  angular  momentum  of  the  underlying  physical  problem  as  well  as  its  energy 
evolution.  A  strictly  positive  energy  dissipation,  in  fact  the  exact  energy  dissipation,  is 
obtained  by  design  during  plastic  steps  while  enforcing  the  plastic  consistency  (i.e.  the 
yield  condition)  on  the  final  stress  appearing  in  the  equations  of  motion.  Exact  energy 
conservation  is  attained,  in  particular,  during  elastic  steps.  The  aforementioned  preser¬ 
vation  of  the  plastic  volume  is  obtained  by  a  new  treatment  of  the  geometric  structure 
behind  the  considered  multiplicative  models  of  finite  strain  plasticity.  Namely,  we  present 
a  new  approximation  of  the  reference  and  elastic  metrics  whose  contractions  with  the  in¬ 
cremental  total  and  elastic  strains  lead  exactly  to  the  increment  of  the  total  and  elastic 
natural  volumetric  strains,  respectively.  The  elastic  metric  is  defined  in  the  intermediate 
configuration,  defined  itself  by  the  proper  discrete  approximation  in  time  of  the  plastic 
deformation  gradient  Fp.  The  new  algorithm  extends  to  the  plastic  range  existing  energy- 
momentum  conserving  schemes  for  nonlinear  elastic  problems,  but  incorporating  a  new 
modified  elastic  stress  formula  consistent  with  this  new  geometric  setting.  The  inherited 
conservation  laws  of  momenta  and,  especially,  the  non-negative  character  of  the  energy 
dissipation  leads  to  an  improved  performance  over  existing,  more  classical  schemes  show¬ 
ing  numerical  instabilities  in  the  considered  highly  nonlinear  geometric  setting  of  large 
deformations  and  strains.  Several  numerical  simulations  are  presented  illustrating  these 
properties. 


8.2.  Appendix  II:  Assumed  Strain  Finite  Element  Methods  for  Conserving 
Time-Stepping  Algorithms 

This  appendix  presents  a  new  assumed  strain  finite  element  formulation  (or  B-bar  method) 
for  the  locking-free  simulation  of  nearly  incompressible  elastic  and  inelastic  solids  in  the 
finite  deformation  dynamic  range  that  also  preserves  the  conservation/dissipation  proper¬ 
ties  of  the  so-called  energy-dissipative  momentum-conserving  (EDMC)  time-stepping  algo¬ 
rithms.  The  general  setting  of  finite  strain  plasticity  is  considered,  including  hyperelastic 
models  as  a  particular  case.  The  main  motivation  of  this  work  is  to  avoid  the  nonlinear 
numerical  instabilities  observed  in  classical  numerical  schemes  with  unbounded  growth  of 
the  energy  (even  in  the  plastic  case)  by  introducing  the  exact  dissipation/conservation  of 
the  energy  in  the  discrete  system  by  design.  The  incorporation  of  the  conservation  laws 
of  linear  and  angular  momenta,  and  the  preservation  of  the  associated  relative  equilibria, 
is  also  obtained.  The  paper  identifies  the  conditions  that  the  linearized  strain  operator 
(or,  simply,  the  B-bar  operator  as  it  is  usually  known)  has  to  satisfy  for  the  preserva- 
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tion  of  these  properties  in  time.  These  conditions  require  the  definition  of  the  assumed 
strain  operator,  originally  developed  by  with  spatial  considerations  only,  accounting  for 
the  temporal  discretization  in  the  definition  of  the  associated  strain  variations.  As  a  re¬ 
sult,  we  arrive  to  a  fully  discrete  system  in  space  and  time  that  shows  exactly  all  these 
conservation/dissipation  laws  of  the  underlying  physical  system,  including  the  exact  plas¬ 
tic  dissipation  of  the  energy,  with  exact  energy  conservation  for  elastic  steps.  Numerical 
simulations  are  presented  to  illustrate  the  performance  of  the  new  formulation. 


APPENDIX  I 


Volume-Preserving  Energy-Dissipative  Momentum- 
Conserving  Time-Stepping  Algorithms  for  Isochoric 
Plastic  Models  of  Multiplicative  Plasticity 


Based  on  the  papers: 

Armero,  F.  &  Zambrana,  C.  [2007]  “Volume-Preserving  Energy- Momentum 
Schemes  for  Isochoric  Multiplicative  Plasticity,”  Computer  Methods  in 
Applied  Mechanics  and  Engineering,  196,  4130-4159. 

Armero,  F.  &  Zambrana,  C.  [2006]  “Numerical  Integration  of  the 
Nonlinear  Dynamics  of  Elastoplastic  Solids,”  Proceedings  of  the  III 
European  Conference  on  Computational  Mechanics  (ECCM-06),  Lis¬ 
bon,  Portugal. 
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1.1.  Introduction 

Classical  time-stepping  algorithms  for  the  numerical  integration  of  problems  in  solid 
and  structural  dynamics,  like  the  Newmark  scheme  Newmark  [1959]  and  its  variations 
(see  Hughes  [1987]  for  a  complete  account),  have  been  shown  to  exhibit  serious  limita¬ 
tions  when  applied  to  geometrically  nonlinear  problems  despite  their  good  performance 
in  the  linear  setting.  These  limitations  include  not  only  the  lack  of  preserving  impor¬ 
tant  conservation  laws  of  the  motion,  like  the  conservation  law  of  angular  momentum, 
but  also  to  numerical  instabilities  in  the  form  of  an  unbounded  growth  of  energy.  These 
instabilities  limit  severely  the  time  step  size  that  can  be  used  in  the  simulations,  even 
in  algorithms  that  have  been  proven  to  be  unconditionally  stable  in  the  linear  range. 
This  situation  has  motivated  the  development  of  the  so-called  energy-momentum  schemes, 
where  the  conservation  laws  of  momenta  and  energy  are  embedded  in  the  algorithm  by 
design;  see  Crisfield  &  Shi  [1994], Gonzalez  [2000], Simo  &  Tarnow  [1992],  among 
others,  and  Armero  &  Romero  [2001a], Armero  &  Romero  [2001b]], Kuhl  &  Cr¬ 
isfield  [1997], Kuhl  &  Ramm  [1996]  for  extensions  exhibiting  a  controllable  numerical 
dissipation  in  the  high-frequency  to  handle  the  high  numerical  stiffness  of  the  systems  of 
interest.  We  note  in  this  respect  that  classical  high-frequency  dissipative  schemes,  like  the 
HHT  (see  HUGHES  [1987]),  loose  also  their  unconditional  stability  in  the  nonlinear  range 
Armero  &  Romero  [2001a]. 

These  numerical  instabilities  have  also  been  observed  in  the  elastoplastic  range  at 
finite  strains,  thus  motivating  the  extensions  of  the  above  energy-momentum  schemes  for 
elastic  systems  to  this  inelastic  setting,  as  considered  in  Meng&  Laursen  [2002], Noels 
ET  AL  [2004].  A  major  challenge  in  these  problems  is  the  proper  integration  of  the  plas¬ 
tic  evolution  equations  for  the  plastic  internal  variables  giving  also  the  stresses  or,  in 
other  words,  the  proper  return  mapping  algorithm.  In  this  context,  the  work  presented  in 
Meng&  Laursen  [2002]  considers  standard  existing  return  mapping  algorithms  with  an 
additional  projection  step  on  the  stress,  along  the  lines  presented  in  Gonzalez  [2000]  for 
elastic  problems,  to  recover  the  proper  energy  dissipation  in  the  plastic  range.  This  step, 
however,  disturbs  the  satisfaction  of  the  initially  enforced  plastic  consistency  condition  for 
the  final  resulting  stress.  The  algorithms  presented  in  Noels  ET  al  [2004]  considered 
hypoelastic  based  models,  thus  being  able  to  impose  only  the  conservation/dissipation  in 
very  specific  particular  cases. 

We  have  recently  presented  in  Armero  [2006], Armero  [2005]  an  alternative  strategy 
leading  to  new  energy-dissipative  momentum-conserving  time-stepping  algorithm  for  mul¬ 
tiplicative  plasticity  ( F  =  F'FP)  that  exhibits  the  exact  energy  dissipation  of  the  underly¬ 
ing  physical  system,  while  enforcing  exactly  the  consistency  condition  on  the  final  stresses. 
The  development  of  the  new  algorithm  follows  the  same  arguments  leading  to  these  con¬ 
servation/dissipation  properties  for  the  continuum  problem,  identifying  in  the  process  the 
proper  discrete  approximation  of  the  plastic  strain  rate  that  preserves  the  dissipative  char¬ 
acter  of  t  he  (discrete)  plastic  flow.  The  final  scheme  is  implemented  following  the  classical 


Final  Report.  FA9550-05-1-0117 


23 


structure  of  return  mapping  algorithms  consisting  of  an  elastic  trial  state  followed  a  plastic 
corrector  imposing  the  yield  condition  on  the  stresses  (or  its  viscoplastic  regularization). 
However,  in  contrast  with  existing  exponential  return  mapping  algorithms,  originally  pro¬ 
posed  in  Eterovich  &  Bathe  [1990],Cuitinho  &  Ortiz  [1992], Simo  [1992], Weber 
&  Anand  [1990]  (see  complete  details  in  the  monograph  SlMO  [1998]),  the  new  scheme 
involves  an  algebraic  approximation  of  the  flow  rule  in  the  updated  plastic  deformation 
gradient  F^+ , ,  not  enforcing  automatically  the  isochoric  character  of  the  plastic  flow  in 
isochoric  plastic  models  like  the  models  of  J2-flow  theory  for  metals. 

It  is  precisely  the  goal  of  this  contribution  to  develop  alternative  algorithms  that  also 
exhibit  this  volume-preserving  property.  As  shown  in  the  developments  below,  the  cru¬ 
cial  observation  in  the  accomplishment  of  this  goal  is  the  proper  approximation  of  the 
geometric  structure  behind  the  considered  multiplicative  models  of  finite  strain  plasticity. 
In  this  way,  the  new  algorithm  considers  explicitly  the  approximation  of  the  metrics  in 
the  reference  and  intermediate  configurations  that  define  the  change  of  volumetric  strains, 
total  and  elastic  respectively,  in  these  models.  The  actual  definition  of  the  intermediate 
configuration  in  the  discrete  setting  of  a  time  step  requires  special  considerations  for  the 
proper  approximation  of  this  geometric  structure.  Altogether,  we  develop  in  this  work  a 
new  volume-preserving  energy- dissipative  momentum- conserving  scheme  that  exactly  pre¬ 
serves  the  isochoric  character  of  the  plastic  flow  in  the  aforementioned  models  of  metal 
plasticity  while  still  preserving  exactly  the  conservation  laws  of  linear  and  angular  mo¬ 
menta,  and  also  recovering  the  exact  energy  dissipation  of  the  physical  system.  This 
means,  in  particular,  that  exact  energy  conservation  is  obtained  for  elastic  steps,  extend¬ 
ing  to  the  elastoplastic  case  existing  energy-momentum  schemes  with  the  aforementioned 
improved  stability  properties. 

An  outline  of  the  rest  of  the  paper  is  as  follows.  Section  1.2  presents  a  summary  of 
the  governing  equations  of  finite  strain  multiplicative  plasticity  in  the  dynamic  range  with 
an  outline  of  the  arguments  leading  to  the  conservation/dissipation  properties  of  interest 
here,  including  the  geometric  structure  behind  these  equations.  Following  the  approach 
advocated  here,  it  is  the  analog  of  these  very  same  arguments  that  leads  to  the  formu¬ 
lation  of  the  new  volume-preserving  energy-dissipative  momentum-conserving  scheme,  as 
developed  in  Section  1.3.  Section  1.4  present  representative  numerical  simulations  illustrat¬ 
ing  the  performance  of  the  new  algorithm  in  comparison  with  existing  schemes.  Finally, 
Section  1.5  includes  some  final  remarks. 


1.2.  Problem  Definition 

We  summarize  in  this  section  the  problem  of  interest  in  this  work.  In  this  way, 
Section  1.2.1  describes  the  governing  equations,  including  the  equations  of  motion  and  the 
relations  defining  the  multiplicative  models  of  finite  strain  plasticity.  Important  in  our 
developments  is  to  fully  characterize  the  dissipative  character  of  the  final  equations,  as 
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L  p=  Fp Fp' 1  =  G e  [  Dp+  Wp] 

FIGURE  1.2.1  Finite  strain  multiplicative  plasticity.  Problem  definition 

done  in  that  section,  and  the  geometric  setting  defining  the  volumetric  strain  changes  in 
their  elastic  and  plastic  parts.  This  is  discussed  in  Section  1.2.2.  We  refer  to  SlMO  [1998] 
for  complete  details  on  the  results  summarized  here. 

1.2.1.  The  governing  equations  and  energy  dissipation  in  multiplicative 
plasticity 

We  are  interested  in  the  motion  < p  :  B  x  [0,  T]  -a-  M'5  of  a  solid  represented  by  B  c  R3 
in  the  time  interval  T,  determined  by  the  weak  equation 

[  Po<p-SipdV  +  I  S:  (FTGRAD(Sip))s  dV  =  Gext(ip ,  Sip)  ,  (1.2.1) 

JB  JB 

for  all  admissible  variations  S(p,  that  is,  Sip  =  0  on  d^B  c  dB  (the  part  of  the  boundary 
where  the  deformation  ip  is  imposed).  Here  we  have  considered  the  deformation  gradient 
F  —  GRAD (<£>),  the  second  Piola-Kirchhoff  stress  tensor  S  =  ST  (symmetric),  the  refer¬ 
ence  density  pa  and  the  solid’s  acceleration  ip  =  d2ip/dt2.  The  right-hand-side  of  (1.2.1) 
corresponds  to  the  virtual  work  of  the  external  loading  denoted  generically  by  Gext. 

We  are  interested  in  models  of  multiplicative  plasticity  defined  by  the  elastoplastic 
decomposition 

F  =  FeFp  ,  (1.2.2) 

defining  locally  the  intermediate  configuration  Opx  by  the  plastic  deformation  gradient  F'p: 
see  Figure  1.2.1.  The  elastic  part  of  the  deformation  gradient  Fe  defines  the  stresses  S 


through  the  hyperelastic  relations 

-i  -  -T  -  d\Ve 

S  =  Fp  SFP  with  S  =  2-^  ,  (1.2.3) 

in  the  intermediate  configuration  Opx  for  the  elastic  potential  We(Ce)  in  terms  of  the 

T 

elastic  right  Cauchy-Green  tensor  Ce  =  Fe  Fe  (by  frame  indifference).  We  note  the 
relation 

We  =  S:\de ,  (1.2.4) 

giving  the  change  of  elastic  strain  energy. 

Assuming  for  simplicity  the  case  of  no  external  loading  (that  is,  Gext  =  0  and  d^B  = 
0),  the  insertion  of  the  velocity  (fi  (now  an  admissible  variation)  in  the  variation  slot  of 
(1.2.1)  leads  to  the  energy  relation 

4  /  \po\M2  dV  +  We(Ce)  +  H(a)  =  -  f  [ S:Dp  +  qa ]  dV  ,  (1.2.5) 

dt  \_JB  2  J  Jb 

' - v - '  " - V - " 

total  energy  H(t)  plastic  dissipation  T> 

identifying  the  plastic  strain  rate 

Dp  :=  sym  Ge'  Lp  =  ^  (fp~T CF1^  -  Ce)  ,  (1.2.6) 

for  Lp  =  FPFP  *  and  Ge  :=  Ce  \  In  (1.2.5)  we  have  introduced  the  hardening  potential 
H(a)  in  terms  of  the  strain-like  internal  variable  a  and  its  conjugate  stress-like  variable 
q  =  — dH/da .  We  assume  for  simplicity  the  case  of  isotropic  hardening  in  terms  of  the 
scalar  equivalent  plastic  strain  a. 

The  role  of  Ge  in  the  definition  of  the  plastic  strain  rate  Dp  in  (1.2.6)  is  to  be  noted, 
defining  the  proper  geometric  setting  in  the  intermediate  configuration  Ovx.  Indeed,  we 
can  identify  this  (symmetric  positive  definite)  tensor  Ge  =  Ce  as  the  metric  in  that 
configuration,  allowing  to  define  the  “two-cova”  tensor  associated  with  Lp  before  taking 
its  symmetric  part.  Similarly,  we  identify  the  metric  G  =  C~l  for  C  =  FTF  in  the 
reference  configuration  Ox]  see  Figure  1.2.1. 

For  the  developments  below,  it  is  important  to  emphasize  that  the  actual  calculation 
of  the  energy  evolution  (1.2.5)  identifies  the  plastic  strain  rate  (1.2.6)  as  it  appears  in  the 
plastic  dissipation  V.  It  is  precisely  the  physical  requirement 

V  >  0  (second  law)  (1-2.7) 

that  motivates  the  use  of  this  plastic  strain  rate  Dp  in  the  formulation  of  plastic  models. 
A  general  elastoplastic  model  can  be  written  as 

Dp  =  7  N(f)§(S,  q]  Ge)  , 

Wp  =  7  MWP(S.  q]  Ge)  , 

a  =  7  n<t>q{S,q]Ge)  , 


(1.2.8) 

(1.2.9) 

(1.2.10) 
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in  the  intermediate  configuration  Opx.  Here,  we  have  introduced  the  plastic  spin  Wp 


Wp  :=  skew 


(1.2.11) 


in  (1.2.9)  so  the  whole  plastic  deformation  gradient  rate  Fp  is  determined  by  the  plastic 
evolution  equations  (I.2.8)-(I.2.9).  This  requires  nine  components  in  three-dimensional 
problems,  given  by  the  six  components  of  the  (symmetric)  plastic  strain  rate  Dp  and  the 
three  components  associated  to  the  skew  plastic  spin  Wp.  Equation  (1.2.8)  corresponds 
to  the  flow  rule  whereas  equation  (1.2.10)  defines  the  so-called  hardening  law. 


The  plastic  multiplier  7  in  equations  (I.2.8)-(I.2.9)  is  determined  by  the  Kuhn- Tucker 
loading/unloading  and  consistency  conditions 


(f>(S,q-,Ge)<  0,  7>0,  fa  =  0  ,  fa  =  0  (1.2.12) 


for  a  general  yield  function  (f)(S,  q;  Ge).  Note  the  dependence  on  the  metric  Ge  of  the 
intermediate  configuration  Opx  of  both  the  yield  function  and  the  plastic  flow  vectors  in 
(I.2.8)-(I.2.9),  besides  the  stress  variables  S  and  q.  The  Perzyna  regularization  of  equations 
(1.2.12)  defines  the  viscoplastic  model 


7 


g(<t>) 


for  g{<j>)  = 


0 


for 

for 


<  0 
>  0 


(1.2.13) 


i)  '  “ ' ' '  [  monotonically  increasing 

and  a  viscous  parameter  77  >  0,  recovering  the  elastoplastic  equations  (1.2.12)  in  the  limit 
rj  — >  0. 


A  typical  choice  for  the  flow  vectors  in  (I.2.8)-(I.2.9)  is  given  by  the  so-called  associated 
relations 

N*s  =  %  and  »  (L2-14) 

with  a  vanishing  plastic  spin  Mwp  =  0,  especially  for  isotropic  models.  The  flow  vectors 
(1.2.14)  result  also  from  the  classical  principle  of  maximum  plastic  dissipation;  see  SlMO 
[1998], 


Remark  1.2.1  The  physical  condition  (1.2.7)  is  easily  satisfied  by  the  associated  plastic 
flow  vectors  (1.2.14)  for  yield  functions  of  the  form 

<HS,  q;G‘)=  J(S;  G«)  -  [<r„  -  ?]  ,  (1.2.15) 

for  the  initial  yield  limit  oy  >  0  and  a  positively  homogeneous  function  of  degree  one 
(i.e.  (j)(XS]  G'')  =  \({)(\S:  Grt)  for  A  >  0).  After  using  Euler’s  theorem  of  homogeneous 
functions,  the  plastic  dissipation  reads  then, 

V  =  J  \j\  7  °y  dV  >0  ,  (1.2.16) 

by  the  Kuhn-Tucker  condition  (1.2.12)4.  A  similar  argument  gives  the  non- negative  char¬ 
acter  of  the  dissipation  in  (1.2.5)  for  the  viscoplastic  case  (1.2.13).  □ 
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Remark  1.2.2  The  equation  of  motion  (1.2.1)  has  also  additional  conservation  laws  of 
momenta  given  the  symmetries  in  the  resulting  dynamical  system.  In  this  way,  for  the 
case  of  no  external  loading  ( Gext  =  0  and  7F B  =  0),  we  have 


l 


Po<p  dV 


'  B 


constant  and  j 


x  pQ(p  dV  =  constant  , 


(1.2.17) 


along  the  solutions  of  the  problem.  The  vectors  in  (1.2.17)  correspond  to  the  linear  and 
angular  momentum,  respectively.  Their  conservation  follows  easily  by  inserting  the  varia¬ 
tions  Sep  =  c  and  6(fi  =  c  x  ip  for  all  c  e  R3  in  (1.2.1),  admissible  variations  for  the  case 
considered  here.  □ 


1.2.2.  The  volume  change  and  isochoric  plastic  models 

Metals  are  known  to  exhibit  no  plastic  change  of  volume  in  the  strain  ranges  of  interest. 
In  the  finite  deformation  setting  considered  in  the  previous  section,  the  natural  volumetric 
strain  is  given  by 

ev:=\ogJ  for  J  =  det  [F]  ,  (1.2.18) 

the  Jacobian  J  and  its  natural  logarithm  log .  The  multiplicative  decomposition  (1.2.2) 
leads  to  the  additive  decomposition 

£v=el+ep,  (1-2.19) 

for  the  elastic  and  plastic  parts 

el  =  log  Je  and  £%  =  log  Jp  ,  (1.2.20) 

in  terms  of  the  elastic  and  plastic  Jacobians  Je  =  det  [Fe]  and  -Jp  =  det  \FP] .  respectively. 
Fundamental  to  the  developments  below  are  the  rate  relations 

£v  =  \d  :G  and  £ev  =  \ce  :  Ge  ,  (1.2.21) 

Zj  Zt 

for  the  total  and  elastic  metrics  G  and  Ge  in  the  reference  Ox  and  intermediate  configu¬ 
rations  Opx ,  respectively.  Given  the  plastic  strain  rate  (1.2.6),  we  can  write 

Ge  :  Dv  =  -C  :  (Fp~XGeFp  T)  -  -Ce  :  Ge  =  ev  -  eev  =  epv  ,  (1.2.22) 

2  V  /  2 

after  noting  the  relation 

G  =  Fp~1GeFp~T  ,  (1.2.23) 

as  a  simple  algebraic  calculation  shows.  The  general  flow  rule  (1.2.8)  leads  then 

%  =  7  N+s  :  , 


(1.2.24) 
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and  hence  we  have 

Jp  =  1  if  N^-  :Ge  =  0  ,  (1.2.25) 

recovering  the  desired  isochoric  character  of  the  plastic  flow  when  this  last  condition  holds. 
It  is  the  goal  of  this  work  to  develop  integration  algorithms  for  the  plastic  evolution  equa¬ 
tions  (I.2.8)-(I.2.10)  that  conforms  with  this  estimate  while  exhibiting  the  energy  dissipa¬ 
tion  (1.2.5)  and  the  momentum  conservation  laws  discussed  in  Remark  1.2.2. 

Relation  (1.2.23)  allows  also  to  define  the  plastic  deformation  gradient  Fp  in  terms  of 
the  metrics  G  and  G  .  Indeed,  introducing  the  polar  decompositions 

F  =  RU  and  Fe  =  ReUe  ,  (1.2.26) 

we  can  write 

Fp  =  Gpl/2A  G  1/2  for  the  rotation  A  =  ReT R  ,  (1.2.27) 

as  a  simple  calculation  shows. 


Remark  1.2.3  A  typical  example  of  isochoric  plastic  model  is  given  by  the  von  Mises 
yield  function,  written  in  the  form  (1.2.15)  with 


0(5;  Ge)  =  \Jgc~x  DEVg,.[S]  :  DEVfJe[S]  G'  1 

=  y/ dev  [r]  :  dev  [r]  =  ||  dev  [r]  ||  ,  (1.2.28) 


where  we  have  introduced  the  deviatoric  part  of  the  stress  S  in  the  intermediate  configu¬ 
ration  Opx 

DEVGe[S]  =S-~  (S  :  Ge_1)  Ge  ,  (1.2.29) 

in  the  metric  Ge  of  that  configuration,  and  the  corresponding  classical  deviatoric  part  in 
the  current  configuration 

dev  [t]  =  Fe  DEVGe[S]  FeT  =  r  -  \  (r  :  1)  1  ,  (1.2.30) 

O 

of  the  Kirchhoff  stresses  r  =  FSFT  =  FeSFe  .  The  associated  plastic  flow  vector  is 
then  given  by 


N, 

*s  d  S 


i  Ge  1  DEVGe[S'l  Ge 

0 


dev  [r] 
1 1  dev  [t] 


(1.2.31) 


which  can  be  easily  seen  to  satisfy  the  condition  (1.2.25). 


□ 


1.3.  Volume-Preserving  Energy-Dissipative  Momentum- Conser¬ 
ving  Schemes 

We  develop  in  this  section  the  new  EDMC  scheme  for  the  integration  of  the  plastic¬ 
ity  problem  summarized  in  the  previous  section.  Section  1.3.1  presents  general  energy- 
dissipative  momentum- conserving  approximations  in  this  context,  while  Section  1.3.2  dis¬ 
cusses  approximations  of  the  corresponding  geometric  structure  that,  in  addition,  preserves 
the  volume  change  of  the  plastic  flow.  Section  1.3.3  discusses  briefly  the  numerical  imple¬ 
mentation  of  the  resulting  scheme. 


1.3.1.  Energy-dissipative  momentum-conserving  approximations 

The  interest  here  is  the  development  of  one-step  time-stepping  algorithms  for  the 
solution  of  the  governing  equation  (1.2.1).  To  this  purpose,  we  consider  the  general  ap¬ 
proximation 

Ipn+1  -  Lpn  =  At  u*  (1.3.1) 

Vn+1^Vn  -^  +  S,:  (if  GRAD  {Sip)) s  dV  =  GexU  (1.3.2) 

for  all  admissible  variations  Sip  in  a  generic  time  step  [tn,tn+ 1]  with  At  =  tn+\  —  tn 
and  a  generic  second-order  approximation  of  the  external  loading  (say  the  original  Gext 
evaluated  at  the  mid-point  (tn+tn+ 1)/2).  The  quantities  (•)*  need  to  be  defines  such  that 
the  conservation/dissipation  laws  described  in  the  previous  section  are  inherited  by  the 
discrete  equations  (I.3.1)-(I.3.2). 

The  conservation  of  linear  momentum  (ln+ 1  =  ln)  for  no  external  loading  is  triv¬ 
ially  satisfied  as  it  corresponds  to  the  discrete  equation  (1.3.1).  Following  an  argument 
completely  analogous  to  the  one  discussed  in  Remark  1.2.2  leading  to  the  conservation 
laws  (1.2.17)  for  the  continue,  we  obtain  the  conservation  law  of  angular  momentum  (i.e. 
jn+i  =  jn  for  no  external  loading)  for  a  symmetric  stress  tensor  S*  and  the  choice 

F*  =  Fn+ 1  =  ^  (Fn  +  Fn+ 1)  and  u*  ||  un+i  =  \  (vn  +  vn+l)  ,  (1.3.3) 

as  an  algebraic  calculation  shows;  see  SlMO  &  Tarnow  [1992],  Armero  &  Romero 
[2001a].  The  fact  that  the  velocity  approximation  u*  needs  only  to  be  a  vector  parallel  to 
the  mid-point  value  vn+i  has  been  exploited  in  Armero  &  Romero  [2001a],  Armero  & 
Romero  [2001b]]  as  a  way  to  introduce  numerical  dissipation  to  handle  the  high  numerical 
stiffness  of  the  system  of  equations  (I.3.1)-(I.3.2). 

Similarly,  following  the  same  arguments  leading  to  the  energy  evolution  equation 
(1.2.5),  we  obtain  for  the  case  of  no  external  loading  ( Gext  =  0  and  dlfB  =  0)  the  in¬ 
cremental  energy  relation 

Hn+1-Hn  =  J^  (W‘+1  -  Wt)  -  S,  :  IaC'  dV 
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’  B 


T( 


-T  _ _ -1 


FP  A CFP  -  A Ce  )+q*Aa 


dV  ,  (1.3.4) 


T> 


for  the  total  energy  H 3  (i  =  {0, 1})  defined  in  (1.2.5),  after  introducing  the  stress 


S*  =  F!>S*Fl>  , 


(1.3.5) 


for  some  approximation  Fl’  defining  the  discrete  intermediate  configuration  OvXif , 
hardening  variable 


TLn+i  -  Hn 

q*  = - 


®n+ 1 


and  the 
(1.3.6) 


so  g*  Aa  =  AH.  Here,  we  have  used  the  notation  A(-)  =  (-)n+i  —  (-)n- 


A  first  step  in  recovering  the  continuum  estimate  (1.2.5)  consists  in  the  proper  approx¬ 
imation  of  the  stress  formula  (1.2. 3)2  so  the  first  integral  of  the  right-hand-side  of  (1.3.4) 
vanishes  (or  it  is  negative  to  account  for  numerical  dissipation;  see  Remark  1.3.1  below). 
Here  we  consider  the  discrete  formula 


S*  =  S(Cen+i)  +  2 


We(C°n+1)  -  W%Cen)  -  S(Cen+i)  :  \AC 


G%ACe 


A  C«G% 


G%ACeGl 


(1.3.7) 


for  an  approximation  of  the  elastic  metric  G%  as  developed  in  the  next  section.  Here 
S(Ce  1)  denotes  the  gradient  formula  (1.2. 3)2  evaluated  at  Cc:'i  =  {C^  +  CA  i)/2.  The 
stress  formula  (1.3.7)  can  be  seen  to  be  a  projection  of  this  value  to  give  the  incremental 
elastic  strain  energy  when  contracted  with  |ACe.  It  is  a  modification  of  the  formula 
originally  proposed  in  GONZALEZ  [2000]  accounting  for  the  presence  of  the  elastic  metric 
G%,  hence  leading  to  a  properly  invariant  formula.  We  note  that  the  denominator  in  (1.3.7) 
only  vanishes  when  A Ce  does  since 

G%ACe  :  A CeG%  =  ||(?f/2AC'eG*1/2||2  >  0  ,  (1.3.8) 


given  the  positive-definite  character  of  the  metric  G%.  Therefore,  the  second-term  in  (1.3.7) 
vanishes  when  A  Ce  =  0. 


The  key  aspect  of  the  discrete  energy  evolution  (1.3.4)  is  that  it  reveals  the  proper 
approximation  of  the  plastic  strain  rate  (1.2.6)  so  the  exact  energy  dissipation  is  recov¬ 
ered.  Indeed,  following  the  same  arguments  as  for  the  continuum  problem  we  consider  the 
discrete  plastic  evolution  equations 


skew 


|  [FT  TA CFP  1  -  ACe] 

Gf‘ 


Oin-\-\  Oin 


Ay  N4s(S*,q*;G%) 

A 7  Mwp (S*.q*]  G%) 
A  7  n<t,q(S*,q*]  G%)  j 


(1.3.9) 
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for  the  discrete  plastic  multiplier  A7  defined  by  the  loading/unloading  conditions 

0*  =  0(5*,  g*;  G%)  <  0  ,  Ay  >  0  Ay  </>*  =  0  .  (1.3.10) 

As  in  the  continuum  case,  these  relations  are  replaced  for  a  viscoplastic  model  of  the  form 
(1.2.13)  by 

Ay  ,  (1.3.11) 

V 

where  is  defined  by  (1.3. 10)  1. 

The  flow  vectors  in  the  left-hand-side  of  the  discrete  plastic  evolution  equations  (1.3.9) 
are  defined  by  the  assumed  plastic  model.  The  inclusion  of  these  equations  in  the  energy 
evolution  equation  (1.3.4)  leads  to  the  exact  energy  dissipation  for  the  discrete  equations. 
In  particular,  for  the  associated  plastic  models  based  on  yield  functions  of  the  form  (1.2.15) 
we  recover 

V  =  J  Ay  ay  dV  >  0  ,  (1.3.12) 

thus  recovering  (1.2.16),  and  resulting  in  the  exact  energy  conservation  for  elastic  steps 
(Ay  =  0).  We  note  that  these  dissipation/conservation  properties  hold  for  any  approx¬ 
imation  of  the  metric  G%  and  the  plastic  deformation  gradient  F*  defining  the  discrete 
intermediate  configuration  Opx*.  A  second-order  approximation  is  given  by  the  values 

K  =  )  (K  +  f„P+i)  and  G l  =  (C'+1)_1  ,  (1.3.13) 

as  originally  proposed  in  Armero  [2006].  We  develop  alternative  definitions  of  these 
quantities  that  lead  to  the  proper  treatment  of  the  plastic  volume  in  the  next  section. 

Remark  1.3.1  We  refer  to  Armero  [2006]  for  a  discussion  on  additional  considerations 
for  the  inclusion  of  a  controllable  high-frequency  energy  dissipation  along  the  lines  of  the 
EDMC  methods  proposed  in  Armero  &  Romero  [2001a],  Armero  &  Romero  [2001b]] 
for  nonlinear  elastodynamics.  □ 


Remark  1.3.2  For  simplicity  in  the  presentation,  we  have  considered  the  semi-discrete 
in  time  equations  (I.3.1)-(I.3.2).  However,  the  fully  discrete  equations,  accounting  for  the 
spatial  discretization,  need  to  be  considered.  The  basic  displacement  finite  element  model 
is  trivially  recovered  from  (1.3. 1)-(I.3.2).  However,  this  basic  approach  is  known  to  lead 
to  locking  in  the  considered  elastoplastic  context,  as  first  noted  in  Nagtegaal  et  al 
[1974].  The  simulations  presented  in  Section  1.4  consider  the  mixed  treatment  presented 
in  Gonzalez  [2000]  adapted  to  the  elastoplastic  problem  considered  here.  Alternative 
assumed  strain  treatments  preserving  the  conservation/dissipation  properties  of  the  mo¬ 
menta  and  energy  can  be  found  in  Armero  [2008].  The  very  same  conservation/dissipation 
properties  of  the  time  discrete  schemes,  the  interest  in  this  paper,  apply  to  these  cases.  □ 
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1.3.2.  The  discrete  geometric  structure 

It  remains  to  define  the  proper  geometric  structure  of  the  discrete  equations  so  the 
volumetric  properties  of  the  continuum  plastic  flow  are  preserved.  The  key  aspect  is  given 
by  the  relations  (1.2.21),  whose  discrete  analog  we  seek.  That  is,  we  require  that  the 
discrete  metrics  G *  and  G%  in  the  reference  and  intermediate  configurations,  respectively, 
satisfy  the  incremental  relations 

Aev  =  ^AG  :  G*  and  Ae%  =  ^ACe  :  G%  ,  (1.3.14) 

for  the  increments  of  total  and  elastic  volumetric  strains  in  (1.2.20),  while  defining  a  second- 
order  approximation  of  the  continuum  values  C_1  and  Ce  ,  respectively. 

Following  the  same  strategy  as  for  the  stresses  in  the  previous  section,  we  consider 
a  mid-point  approximation  of  these  values  with  a  projection  term  enforcing  the  relations 

(1.3.14).  In  this  way,  we  define 

(1.3.15) 

for  the  metric  in  the  reference  configuration  Ox *,  with  Cn+ 1  =  ( Cn  +  Cn+ 1)/2  and  its 
inverse  C~]_  x ,  Jn+\  =  (det  [Cn+ 1})1^2  and  Jn  =  (det  [Cn])1^2  and 

n  i  2 

(1.3.16) 

for  the  metric  in  the  intermediate  configuration  Opx ,  with  G^+1  =  (C(a  +  Cen+ 1  )/2  and  its 

inverse  C^+l ,  J^+i  —  (det  [C7®+1]) '  '''2  and  .J'n  =  (det  [C);])1//2.  We  note  again  that  the  de¬ 
nominators  in  these  expressions  are  well-defined,  vanishing  when  the  respective  numerators 
vanish,  that  is,  when  AC  and  A Ce  vanish,  respectively. 

The  use  of  the  metrics  (1.3.15)  and  (1.3.16)  in  combination  with  the  discrete  flow  rule 
(1.3.9)  i  leads  directly  to  the  relation 

Ay  N+s  :  G%  =  ^AC  :  ( Ff'GIFf T)  -  ±ACe  :  G% 

=  \  aC  :  G*  -  \  ACe  :  G%  =  As  -  Aee  =  Aep  ,  (1.3.17) 

Li  At 

the  complete  analog  of  the  continuum  relation  (1.2.24),  as  long  as  (1.2.23)  holds  in  this 

^  -i  -T 

discrete  setting,  that  is,  as  long  as  G*  =  Fl'  G%Fl'  .  This  is  obtained  by  considering 
the  approximation  of  the  continuum  relation  (1.2.27)  given  by 

(1.3.18) 


Final  Report.  FA9550-05-1-0117 


33 


The  rotation  A*  corresponds  to  a  second-order  interpolation  of  the  rotation  A  in  (1.2.27) 

T  T 

from  the  end  values  An  =  Ren  Rn  and  An+ 1  =  R(f+lRn+1.  Formula  (L3.18)2  can  be 
efficiently  implemented  (with  algebraic  operations  only)  using  quaternions: 


Qn  A  Qn-\- 1 
Qn  A  Qn+1 


(1.3.19) 


for  the  quaternions  q  associated  to  the  different  rotations  A;  see  Shoemake  [1985].  We 
refer  to  Armero  &  Zambrana  [2007]  for  complete  details,  including  the  linearization  of 
the  rotation  interpolation  formulas  (1.3.18)  and  (1.3.19). 

These  relations  define  completely  the  new  volume-preserving  energy-dissipative  momentum- 
conserving  scheme  and  prove  the  following  theorem. 


Proposition  1.3.1  The  one-step  time-stepping  algorithm  defined  by  the  global  relations 
(I.3.1)-(I.3.2)  with  the  stresses  S *  defined  by  (1.3.5)  in  terms  of  S*  given  by  the  elastic 
stress  formula  (1.3.7)  and  the  approximation  (I.3.9)-(I.3.10)  (or  (I.3.9)-(I.3.11)  for  the 
viscoplastic  case)  with  the  hardening  variable  defined  by  (1.3.6)  satisfy: 

i.  The  conservation  laws  of  linear  and  angular  momentum  are  exactly  preserved. 

ii.  The  energy  evolves  following  the  relation  (Gext  =  0) 

Hn+1  -  Hn  =  -V  ,  (1.3.20) 

for  the  exact  physical  dissipation  V  >  0.  with.  V  =  0  for  elastic  steps. 

Hi.  The  plastic  volume  is  exactly  preserved,  that  is, 

Jp+1  =  Jp  =  1  for  isochoric  plastic  models,  ( 1.3.21 ) 


with  the  last  property  holding  for  the  approximations  (1.3.15)  and  (1.3.16)  of  the  reference 
and  elastic  metrics,  respectively,  the  latter  defined  in  the  intermediate  configuration  Opx # 
determined  by  the  plastic  deformation  gradient  (1.3.18). 


Remark  1.3.3  Other  formulas  for  the  discrete  metrics  (1.3.15)  and  (1.3.16)  satisfying  the 
relations  (1.3.14)  are  possible.  For  example,  the  projection  terms  can  be  defined  in  terms  of 
the  metrics  themselves,  giving  an  implicit  definition  of  these  metrics.  Similarly,  the  stress 
projection  term  in  the  elastic  stress  formula  (1.3.7)  can  be  defined  with  the  value 
instead  of  the  metric  G(  and  still  results  in  all  the  properties  summarized  in  Theorem 
1.3.1.  □ 
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Remark  1.3.4  We  note  that  the  volume-preserving  property  (1.3.17)  is  not  only  impor¬ 
tant  for  the  isochoric  plastic  models  as  emphasized  in  (1.3.21),  but  also  in  models  where  the 
plastic  volume  change  is  an  important  physical  aspect  to  capture  correctly,  like  in  dilatant 
models  of  soil  mechanics.  □ 


1.3.3.  The  numerical  implementation 

The  numerical  solution  of  the  discrete  plastic  evolution  equations  (1.3.9)  can  be  ap¬ 
proached  through  the  common  structure  of  return  mapping  algorithms  consisting  of  an 
elastic  trial  state  followed  by  a  plastic  corrector  if  necessary.  Some  calculations  show  that 
the  solution  for  an  elastic  step  Ay  =  0  is  given  by  F f  =  F£.  The  plastic  corrector  re¬ 
quires  the  solution  of  the  system  of  equations  (1.3.9)  for  {F^+1,an+i}  while  evaluating  the 
stresses  S*  through  (1.3.7)  and  the  hardening  variable  q*  through  (1.3.6). 

This  solution  is  efficiently  accomplished  through  a  two-level  scheme  with  two  nested 
Newton  iterations.  The  “upper  level”  iteration  drives  the  residual 

■hy>(^n+i)  =  (oin+i ),  q* (ura+i ) i  G * (crn-|-i ))  (1.3.22) 


to  zero,  solving  for  the  updated  equivalent  strain  an+ 1  while  enforcing  the  consistency 
condition  (or  the  equivalent  relation  for  the  viscoplastic  problem).  The  evaluation  of  the 
stresses  and  internal  variables  for  a  fixed  value  of  a  involves  the  evaluation  of  the  plastic 
flow  and  plastic  spin  rules  (1. 3. 9)1,2-  These  equations  are  solved  for  the  updated  value 
F%+ 1  through  the  “lower  level”  Newton  iteration  driving  the  residual 


RMK+i) 


Ff\cn+ 1  -  Cn)Ff '  -  (C«+1  -  Cl)  -  2A7JV^ 


AXIAL 


G:  '  (F; 


<  P 

n+1 


Ay  M\\rp 


(1.3.23) 


to  zero  with  Ay  =  (o:n+i  —  an )  Afo,  for  the  fixed  value  an+ 1  at  this  lower  level.  Here  we 
used  the  notion  of  the  axial  vector  of  a  tensor  through  its  skew  part  (denoted  by  axial  [•]), 
so  the  residual  in  (1.3.23)  has  nine  components  for  the  nine  components  of  F%+1  after 
noting  the  symmetry  of  the  first  component  of  the  residual  in  this  equation.  We  refer 
to  Armero  &  Zambrana  [2007]  for  complete  details  of  the  numerical  implementation, 
including  the  consistent  tangents  in  the  solution  of  the  problems  (1.3.22)  and  (1.3.23),  as 
well  as  the  global  algorithmic  consistent  tangent  used  in  the  construction  of  the  stiffness 
matrix  in  typical  Newton-Raphson  solutions  of  the  global  equations  (I.3.1)-(I.3.2). 


FIGURE  1.3.1  Tumbling  L-shaped  block.  Initial  configuration  and  loading  with  finite  element 
mesh  (117  constant  volume  mixed  trilinear  bricks). 


1.4.  Representative  Numerical  Simulations 

To  illustrate  the  numerical  properties  of  the  newly  developed  volume-preserving  energy- 
dissipative  momentum-conserving  scheme  (referred  simply  in  what  follows  as  the  EDMC- 
VP  scheme)  we  consider  the  problem  of  a  tumbling  L-shaped  block,  first  considered  in 
SiMO  &  Tarnow  [1992]  in  the  elastic  range  and  in  Meng&  Laursen  [2002], Noels  et 
AL  [2004]  in  the  context  of  elastoplasticity.  Figure  1.3.1  depicts  the  initial  configuration  of 
the  block,  with  its  geometry  and  loading.  This  consists  of  point  loads  at  the  nodes  of  the 
block  bases  given  by 

(  Pc  x  t  for  0  <  t  <  25, 

pit)  =  \  Po  x  (50  -  t )  for  25  <  t  <  50,  (1.4.1) 

\  0  for  t  >  50, 

for  the  vector  pa  =  [4  10  12] T  and  the  opposite  vector  in  the  opposite  base  of  the  block,  in 
the  Cartesian  system  defined  by  the  three  orthogonal  directions  along  the  block  edges.  The 
block  is  then  in  free-flight  after  an  initial  period  of  loading.  Figure  1.3.1  depicts  also  the 
assumed  finite  element  discretization,  consisting  of  117  brick  elements.  The  mixed  finite 
element  strategy  indicated  in  Remark  1.3.2  has  been  considered  with  a  piece-wise  constant 
interpolation  for  the  volume  in  combination  with  a  trilinear  displacement  interpolation. 

We  consider  J2-flow  theory  for  the  material  response,  characterized  by  the  von  Mises 
yield  function  (1.2.28)  and  the  associated  flow  vectors  (1.2.14).  Linear  isotropic  hardening  is 
assumed,  that  is,  y(a)  =  ay  —  q  =  ay  +  Ha  for  the  initial  yield  limit  ay  and  linear  hardening 
modulus  H.  We  consider  Hencky’s  law  for  the  hyperelastic  response,  characterized  by  the 
elastic  potential 

1  3 

We(Ce)  =  -K(logJe)2  +/(  ^  (£a)2  , 

A=  1 


(1,4.2) 
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FIGURE  1.3.2  Tumbling  L-shaped  block.  Deformed  configurations  with  the  distribution  of 
the  equivalent  plastic  strain  a  during  the  early  stages  of  the  motion  computed  with  the  new 
EDMC-VP  scheme. 
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FIGURE  1.3.3  Tumbling  L-shaped  block.  Temporal  evolution  of  the  three  components  of  the 
angular  momentum  (left)  and  the  kinetic,  potential  and  total  energies  (right)  for  the  solution 
computed  with  the  new  EDMC-VP  scheme. 


for  the  bulk  and  shear  moduli  k  and  //.  respectively,  and  the  regularized  logarithmic  prin¬ 
cipal  strains 


£a  ••=  log 


4  =  1,2, 3, 


(1-4.3) 


where  \eA  are  the  elastic  principal  stretches  defined  by  the  eigenvalues  of  Ce  (i.e.  Cr N\  = 
(A'jJ  Nu 4,  no  sum  in  A).  Table  1.4.1  includes  the  values  of  the  different  material  parameters 
assumed  in  the  simulations  presented  here. 


We  consider  several  runs  at  constant  time  steps  At.  Figure  1.3.2  includes  several 
snapshots  of  the  solution  computed  with  the  new  ENDC-YP  scheme  with  At  =  0.5  in  the 
early  stages  of  the  motion.  The  distribution  of  the  equivalent  plastic  strain  a  is  depicted 
over  the  deformed  configuration,  thus  showing  the  extend  of  the  plasticity.  The  large 
displacements  and  strains  involved  in  the  solution  are  to  be  noted. 

Figure  1.3.3  shows  the  evolution  of  the  three  components  of  the  angular  momentum 
j  and  the  energy,  including  the  kinetic  and  potential  energies,  for  the  solution  depicted  in 
Figure  1.3.2,  that  is,  the  solution  computed  with  the  new  EDMC-VP  scheme  and  a  constant 
time  step  of  At  =  0.5.  The  potential  energy  consists  of  the  elastic  strain  energy  plus 
the  hardening  potential  integrated  over  the  solid.  The  exact  conservation  of  the  angular 
momentum  after  the  initial  loading  period  is  confirmed,  as  it  is  the  conservation  of  the 
linear  momentum  (not  shown).  Similarly,  the  monotonic  non- negative  energy  dissipation 
after  this  loading  period  is  confirmed  by  these  plots.  The  plastic  volume  (not  shown)  is 
also  confirmed  to  be  exactly  preserved  at  all  quadrature  points:  Jp  =  1  everywhere  and 
all  times  for  the  assumed  isochoric  model  of  J2-flow  theory. 
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FIGURE  1.4.1  Tumbling  L-shaped  block.  Evolution  of  the  total  energy  in  time  for  the  solu¬ 
tions  computed  with  the  new  EDMC-VP  scheme  (left)  and  the  classical  trapezoidal  rule  (right) 
for  different  time  steps.  The  instabilities  of  the  trapezoidal  rule  is  to  be  contrasted  with  the 
strictly  non-negative  energy  dissipation  of  the  new  EDMC-VP  scheme. 


TABLE  1.4.1  Tumbling  of  a  L-shaped  block.  Material  parameters 


Bulk  modulus 

K 

2500 

Shear  modulus 

P 

500 

Initial  uniaxial  yield  limit 

Oy 

100.0 

Linear  hardening  modulus 

H 

200 

Reference  density 

Po 

100 

Table  1.4.2  shows  the  norm  of  the  residual  during  a  typical  time  step  of  the  global 
Newton-Raphson  scheme  employed  in  the  solution  of  the  global  nonlinear  equations.  An 
asymptotic  quadratic  rate  of  convergence  is  observed,  confirming  the  use  of  the  algorithmic 
consistent  tangent  (see  Armero  &  ZAMBRANA  [2007]).  Table  1.4.3  includes  the  residual 
during  the  nested  Newton  iterations  of  the  return  mapping  in  a  typical  quadrature  point 
undergoing  plastic  deformation.  The  “upper  level”  iteration  enforces  the  consistency  con¬ 
dition,  requiring  the  evaluation  of  the  stresses  with  the  corresponding  update  F%+1  for  a 
given  equivalent  plastic  strain  a.  This  is  done  in  the  “lower  level”  enforcing  the  plastic 
flow  and  spin  rules.  A  quadratic  rate  of  convergence  can  be  observed,  with  a  maximum 
of  four  iterations  per  level,  as  all  the  steps  are  linearized  consistently,  easing  on  the  added 
computational  cost  involved  with  the  new  EDMC-VP  scheme.  In  this  respect,  we  note  that 
a  substantial  part  of  this  added  cost  can  be  traced  to  the  non-symmetry  of  the  algorithmic 
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TABLE  1.4.2  Tumbling  of  a  L-shaped  block.  Convergence  of  the 
Newton-Raphson  scheme  in  a  typical  time  step 


1  2.1521743  •  10+O3 

2  1.6642863  ■  10+O2 

3  1.0519854  •  10+°° 

4  1.9821388  ■  10“O4 

5  5.5203698  ■  lO”08 

6  3.3582068  ■  10”11 


TABLE  1.4.3  Tumbling  of  a  L-shaped  block.  Convergence  of  the  two 
nested  Newton  schemes  to  solve  the  local  return  mapping  equations  in  a 
typical  quadrature  point  for  a  plastic  step.  The  “upper  level”  iteration 
imposes  consistency,  calling  the  “lower  level”  to  evaluate  the  stresses 
and  internal  variables  by  enforcing  the  flow  and  plastic  spin  rules. 


1 

0.33562  •  10+01 

1 

0.12540-  10“01 

2 

0.72429  •  10“04 

3 

0.29601  •  10“08 

4 

0.14922  •  10-12 

2 

0.35841  •  10+°° 

1 

0.14982-  nr02 

2 

0.10431  •  10"05 

3 

0.40102  •  10"11 

3 

0.86176  •  10-°4 

1 

0.36027  •  10"06 

2 

0.81996-  10"13 

4 

0.66080  •  HT11 

1 

0.27060  •  10"13 

consistent  tangent  characteristic  of  energy-momentum  schemes,  even  in  elastic  problems. 
An  efficient  strategy  to  avoid  this  unsymmetry,  relying  on  staggered  symmetric  iterations, 
can  be  found  in  Armero  &  Romero  [2001b]]. 

We  present  in  Figure  1.4.1  a  zoom  of  the  evolution  of  the  total  energy  computed  with 
different  time  steps  (At  =  0.5,  0.75  and  1.00).  In  all  cases,  we  observe  the  monotonic 
dissipation  of  the  new  EDMC-VP  scheme,  confirming  the  exact  energy  conservation  dur¬ 
ing  elastic  steps.  These  solutions  are  to  be  compared  with  the  solution  computed  with 
the  classical  trapezoidal  rule  (Newmark  ryNW  =  1/2  ft  N w  =  1/4)  in  combination  with  an 
exponential  return  mapping  for  the  integration  of  the  plastic  evolution  equations  as  also 


F.  Armero 


40 


shown  in  Figure  1.4.1;  see  SlMO  [1992], SlMO  [1998]  for  a  complete  description  of  the  im¬ 
plementation  of  this  classical  algorithm.  The  instabilities  of  this  scheme  in  the  nonlinear 
range  under  consideration  are  clear:  a  blow  up  of  the  energy  is  observed  at  a  finite  time, 
forcing  the  stoppage  of  the  of  the  simulation.  The  superior,  more  stable,  performance  of 
the  new  EDMC-VP  scheme  is  concluded. 


1.5.  Concluding  Remarks 

We  have  presented  in  this  paper  a  new  energy-dissipative  momentum-conserving  time¬ 
stepping  algorithm  for  finite  strain  multiplicative  plasticity  that  also  preserves  the  volume 
change  of  the  plastic  flow.  In  particular,  it  reproduces  exactly  the  isochoric  plastic  response 
of  isochoric  plastic  models  for  metals.  The  new  scheme  relies  on  a  new  integration  scheme 
of  the  plastic  evolution  equations  (flow  and  plastic  spin  rules  and  hardening  law),  whose 
implementation  can  also  be  recast  in  the  classical  structure  of  return  mapping  algorithms 
consisting  of  an  elastic  trial  step  followed  by  a  plastic  corrector.  The  scheme  is  constructed 
by  following  the  discrete  analogs  of  the  arguments  that  build  the  energy-momentum  dis¬ 
sipation/conservation  properties  in  the  underlying  continuum  problem.  Similarly,  the 
volume-preserving  character  is  constructed  by  the  proper  approximation  of  the  geometric 
structure  behind  the  multiplicative  plasticity  models  of  interest  in  this  work.  The  new 
scheme  thus  developed  exhibits  the  aforementioned  dissipation/conservation  properties  by 
design,  rigorously  and  independently  of  the  time  step  and  model.  In  fact,  the  scheme  is 
completely  general,  applying  to  isotropic  or  anisotropic  elastic  and/or  plastic  models  of  the 
associated  or  non-associated  type.  The  improved  stability  properties  of  the  new  scheme 
have  been  illustrated  for  a  representative  numerical  simulation,  where  standard  schemes 
like  the  trapezoidal  rule  exhibit  numerical  instabilities  in  the  form  of  an  unbounded  growth 
of  the  energy  even  in  the  dissipative  plastic  range  considered  here. 
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II.  1.  Introduction 

Classical  time-stepping  algorithms,  like  the  Newmark  or  HHT  schemes,  developed 
originally  in  the  context  of  linear  elastodynamics,  are  known  to  lead  to  severe  numerical 
instabilities  in  the  nonlinear  finite  deformation  range,  even  for  schemes  that  are  uncondi¬ 
tionally  stable  in  the  linear  range;  see  e.g.  Simo  &  Tarnow  [1992],  Armero  &  Romero 
[2001a],  among  others,  or  the  results  presented  here.  These  instabilities  are  characterized 
by  an  unbounded  growth  of  the  energy,  and  have  been  observed  even  in  the  context  of 
elastoplastic  models  (Meng  &  Laursen  [2002],  Armero  [2006]).  This  situation,  and  the 
lack  of  the  conservation  law  of  angular  momentum  in  many  of  these  classical  schemes,  has 
motivated  a  large  amount  of  recent  literature  on  the  formulation  of  the  so-called  energy- 
momentum  schemes.  These  schemes  inherit  the  conservation  laws  of  energy  and  momenta 
of  the  underlying  physical  system  by  design. 

We  refer  to  Simo  &  Tarnow  [1992],  Crisfield  &  Shi  [1994],  Gonzalez  [2000]  for 
an  illustration  of  energy- momentum  methods  in  nonlinear  elastic  problems,  and  to  Meng 
&  Laursen  [2002],  Noels  et  al  [2004],  Armero  [2006],  Armero  &  Zambrana  [2007] 
for  formulations  considering  the  elastoplastic  range  where  the  goal  is  to  capture  the  exact 
plastic  dissipation  (with  exact  conservation  for  elastic  steps)  while  still  preserving  the 
momentum  conservation  laws.  Extensions  of  these  methods  to  incorporate  an  additional 
controllable  numerical  energy  dissipation  in  the  high-frequency  range  in  order  to  handle  the 
characteristic  numerical  stiffness  of  typical  mechanical  and  structural  system  of  interest 
have  been  proposed  in  Armero  &  Romero  [2001a],  Armero  &  Romero  [2001b]  for 
nonlinear  continuum  elastodynamics  and  in  Armero  &  Romero  [2003],  Romero  & 
Armero  [2002]  in  the  context  of  rod  and  shell  Cosserat  models  of  nonlinear  structural 
dynamics.  We  refer  to  these  time-stepping  algorithms  as  EDMC  schemes  (for  Energy- 
Dissipative  Moment um- Conserving).  They  include,  as  a  particular  case,  some  of  the 
aforementioned  energy- momentum  schemes. 

All  these  references  consider  the  finite  element  method  for  the  spatial  discretization. 
The  consideration  of  a  nearly  incompressible  material  response,  like  the  one  observed 
in  plastically  deforming  metals  and  captured  by  classical  elastoplastic  models  of  J2-flow 
theory,  requires  the  consideration  of  finite  element  formulations  more  sophisticated  than 
the  basic  displacement  model  to  avoid  the  characteristic  volumetric  locking  of  this  basic 
formulation.  To  this  end,  the  so-called  assumed  strain  approach,  as  developed  originally 
by  Nagtegaal  et  al  [1974],  Hughes  [1980]  in  the  infinitesimal  range  and  Simo  et  al 
[1988], Armero  [2000]  in  the  finite  deformation  range  (both  for  static  problems),  becomes 
very  convenient  as  it  only  requires  the  proper  definition  of  the  numerical  approximation  of 
the  strains  and  their  variations  regardless  of  the  material  model  under  consideration.  In 
the  continuum  nearly  incompressible  context  of  interest  here,  these  formulations  are  also 
known  as  “B-bar”  methods. 

B-bar  methods  that  lead  to  locking-free  finite  elements  in  general  configurations  are 
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well-known  by  now,  including  general  linear  and  quadratic  quadrilateral  and  triangular 
elements  for  two  dimensional  problems,  and  similarly  in  3D.  Unfortunately,  their  direct 
consideration  in  the  dynamic  range  of  interest  here  destroys  completely  the  conserva¬ 
tion/dissipation  properties  outlined  above  when  used  in  combination  of  the  aforementioned 
energy-momentum  or  EDMC  schemes.  These  time-stepping  algorithms  rely  on  specific  in¬ 
cremental  properties  of  the  linearized  strain  operator  appearing  in  the  equation  of  motion 
for  a  typical  time  step,  properties  that  a  straightforward  evaluation  of  the  B-bar  operator, 
say  at  the  mid-point,  does  not  have.  This  situation  can  be  traced  back  to  the  nonlinear 
definition  of  the  assumed  deformation  gradient  defining  the  assumed  strain.  For  the  in¬ 
compressible  limit  of  interest  here,  the  assumed  deformation  gradient  involves  a  nonlinear 
scaling  with  its  determinant  or  Jacobian  (another  nonlinear  operation)  and  the  assumed 
Jacobian  defined  through  a  weighted  average  over  the  element. 

All  these  considerations  lead  to  the  need  of  a  new  B-bar  operator  if  the  fundamental 
conservation  laws  of  energy  and  momentum  are  to  be  preserved.  The  new  operator  needs 
to  account  not  only  for  the  discrete  finite  element  interpolations  in  space,  but  also  the 
discrete  structure  in  time  of  the  EDMC  time-stepping  algorithms,  as  presented  in  this 
paper. 


II. 2.  The  governing  equations  and  their  conservation  laws 

We  consider  a  solid  B  C  Rndim  (nciim  =  1,2  or  3)  and  its  motions  ep(X,  t )  in  time  t  for 
the  material  particles  X  E  B,  which  satisfy  the  weak  equation 

I  p0<p-5epdV  +  [  S  :  (FtGRAD[J^])s  dV  =  I  Pob-SepdV  +  [  T  -Sep  dA  ,  (II.2.1) 
JB  JB  '■ - v - '  JB  Jb 

=:±6C{(p,6(p) 

for  all  admissible  variations  Sep,  that  is,  Sep  =  0  on  the  part  of  boundary  d^B  with  imposed 
deformation  ep  =  ep,  complementary  to  the  part  of  the  boundary  ()tB  in  (II. 2.1)  where 
the  tractions  T  are  imposed.  We  have  introduced  in  (II. 2.1)  the  reference  density  of  the 
solid  p0 ,  the  acceleration  ep  =  d2ep/dt 2  ,  the  specific  body  force  b  and  the  second  Piola- 
Kirchhoff  stress  tensor  S,  a  symmetric  tensor  in  the  reference  configuration  B  of  the  solid. 
We  observe  the  appearance  of  the  conjugate  variations  SC/ 2  of  the  right  Cauchy-Green 
tensor  C  =  FrF  for  the  deformation  gradient  F  =  GRAD  [ep]. 

The  particular  form  of  the  governing  equation  (II. 2.1)  leads  to  a  number  of  physical 
conservation  laws,  very  characteristic  of  the  motions  of  solids  and  structures.  In  particular, 
denoting  the  velocity  field  by  V  =  ep,  we  easily  obtain  for  the  case  of  a  free  solid  for  brevity 
(i.e.  for  b  =  0,  T  =  0  and  dvB  =  0)  the  conservation  laws 

l  :=  /  p0V  dV  =  constant  and  j  :=  /  ep  x  pQV  dV  =  constant,  (II. 2. 2) 
Jb  Jb 
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corresponding  to  the  linear  and  angular  momentum,  respectively,  after  using  the  crucial 
properties 

SC(p,c)  =  0  and  SC(p,  c  x  p)  =  0  ,  (II. 2. 3) 

for  all  constant  vector  fields  c  e  M”"liT" .  The  relations  (II. 2. 3)  correspond  to  the  infinitesimal 
generators  of  the  action  of  the  Euclidean  group  ®d''di'"  x  SO(n(\\m)  associated  with  the 
symmetry  of  the  governing  equations  (II. 2.1)  under  translations  and  rotations,  respectively; 
see  e.g.  Marsden  [1992].  This  leads  to  the  existence  of  special  (dynamic)  equilibrium 
solutions  given  by  the  group  motion 


<Pet(X,t)  =  exp  \t  SPiN[J?e]|  pe(X)  +  (J  exp  | rj  spin  [f2e]  |  dr])ve 


(II. 2. 4) 


in  terms  of  two  fixed  vectors  fle  and  Ve,  the  angular  and  translational  velocities  ,  re¬ 
spectively,  and  the  relative  equilibrium  configuration  pe(X)  satisfying  the  equilibrium 
equation 


I  PoOe  x  [f}e  x  pe  T  Ve ]  •  Sp  dV  +  [  S(pe)  :  Fj GRAD  [Sep]  dV  =  0  ,  (II.2.5) 

JB  JB 


again  for  all  admissible  variations  8cp‘,  see  SlMO  ET  AL  [1991].  The  existence  of  these 
relative  equilibria  relies  again  on  the  critical  property  of  the  strain  variations 


SC (cpet,  Sep)  =  SC  (Ve,  EXP  -  t  SPIN  [f2e]  Sp^j 


(II. 2. 6) 


along  the  group  motion  pet(X,t)  in  (II. 2. 4).  Here,  SPIN  [17e]  denotes  the  skew  tensor  with 
axial  vector  !?e,  and  EXP  [spin  [J?e]]  the  rotation  defined  by  the  exponential  map  between 
skew  and  rotation  tensors. 


Finally,  we  note  that  we  always  have  the  relation 


d_ 

dt 


+  W 


dV 


=  :H 


>B 


V  dV  for  V  =  S:-C-W  , 

2 


(II. 2. 7) 


for  a  general  function  W.  Crucial  again  for  obtaining  the  stress  power  in  (II. 2. 7)  is  the 
relation 

SC(p.,p)  =  C,  (II. 2. 8) 

for  the  strain  variations.  The  interest  here  is  the  consideration  of  material  models  with 
W  corresponding  to  the  (internal)  stored  energy  (H  being  the  total  energy  with  the  ki¬ 
netic  energy)  and  D  the  energy  dissipation,  a  non-negative  quantity  by  the  second  law  of 
t  her  mo  dynamics . 

In  particular,  we  are  interested  in  the  case  of  finite  strain  plasticity  characterized  by 
the  multiplicative  decomposition  F  =  FeFp  of  the  deformation  gradient  in  an  elastic  and 
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plastic  part;  see  Armero  [2006]  and  references  therein.  The  stresses  are  then  given  in 


T 

terms  of  an  elastic  potential  IUe(C'e)  for  Ce  =  Fe  Fe  as 

t  dWe 

S  =  FPSFP  =  2  , 

dCe 

(II. 2. 9) 

with  the  plastic  part  Fp  defined  by  the  plastic  evolution  equations  (Lp  = 

pppp~x\ 

Dp  :=  syrn  [CeLp]  =  Fp  TCFp ^  -  Ce  =  7  N^S,  q)  , 

(II. 2. 10) 

Wp  :=  skew  [ CeLp }  =  7  MWP(S,  q)  , 

(II. 2. 11) 

a  =  7  n<j>(S,q)  , 

(II. 2. 12) 

0(<SI,  q)  <  0  ,  7  >  0  ,  7 <f>  =  0and7</>  =  0  , 

(II. 2. 13) 

for  the  yield  surface  cf)(S ,  q)  characterizing  the  elastic  domain  in  stress  space.  We  have 
considered  isotropic  hardening,  modeled  by  the  equivalent  plastic  strain  a  and  the  conju¬ 
gate  stress-like  variable  q  :=  —ddi/da  for  a  hardening  potential  'H(a).  In  this  setting,  the 
internal  energy  W  =  We  +  'H  with  the  plastic  dissipation  given  by  V  =  S  :  Dp  +  qa.  The 

hyperelastic  case  is  recovered  for  a  fixed  Fp.  which  is  the  case  assumed  for  the  relative 
equilibria  (II. 2. 5)  (i.e.  vanishing  of  the  plastic  evolution  or  7  =  0  in  (II.2.10)-(II.2.13)). 


II. 3.  EDMC  time-stepping  algorithms 

The  numerical  solution  of  the  governing  equations  (II. 2.1)  proceeds  with  the  consid¬ 
eration  of  their  spatial  and  temporal  discretizations.  The  spatial  discretization  of  interest 
here  starts  with  the  finite  element  discretizations  of  the  deformation  and  velocity  fields 

ode. 

V(X,tn+i)*Vhn+i(X)=  y  Na(X)  (XA+dA+i)  ,  (II.3.1) 

^  V 

A—l  - v - 

_ nn  A 

—  •^n  +  i 

and 

B  ri  o  d  e. 

V(X,tn+i)  «  Vi+i(X)  =  Y.  NA(X)  vA+,  ,  (II.3.2) 

A=1 

with  i  =  0,  or  1,  and  for  the  shape  functions  NA(X)  associated  to  the  nnode  nodes  with 
nodal  (reference)  coordinates  XA,  displacements  dri  and  velocities  vn+z  in  a  typical  time 
step  [tn,tn+ 1]  with  At  =  tn+ 1  —  tn,  not  necessarily  constant. 

Using  standard  procedures,  together  with  a  one-step  mid-point  interpolation  in  time 
(vn+ 1  :=  ( vn  +  un+i)/2),  we  obtain  the  discrete  algebraic  system  of  equations 

y-M  (vn+1  -vn)+  [  dV  =  fext  ,  (II. 3. 3) 

^ t  Jt3h 

dn+i  dn  At  vn_^_i  , 


(II. 3. 4) 
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for  a  (symmetric)  stress  approximation  S'*  to  be  defined.  Here,  we  have  introduced  the 
standard  nodal  external  force  fext  and  the  mass  matrix  M  defined  as  usual  by  the  assembly 

AHelem  /t  D 

KI'  1  with 

e=l 

MAB  =  [  PoNANB  dV  ,  or  Mab  =  [  p0NAdVSAB  (no  sum)  ,  (II.3.5) 

JBh  J Bh 

e  e 

for  the  consistent  or  lump  forms  of  the  mass.  We  have  also  considered  an  assumed  linearized 
strain  operator  B,  (or,  simply,  the  B-bar  operator)  defined  by  the  relations 

-8C  «  -SC  «  BJd  ,  (II. 3. 6) 

2  2 

for  the  nodal  variations  8d.  Note  the  approximate  signs  in  this  equation,  indicating  nu¬ 
merically  consistent  approximations  (in  fact,  second  order  approximations  of  the  mid-point 
values).  In  particular,  the  stresses  are  assumed  to  be  given  in  terms  of  the  assumed  strains 

_  _ rn _  _ 

C  =  F  F  for  an  assumed  deformation  gradient  F  fa  GRAD[<^]  as  considered  in  the 
following  section. 

The  goal  is  the  development  of  numerical  approximations  that  preserve  the  conservation/dis¬ 
sipation  laws  of  energy  and  momenta  identified  in  the  previous  section  for  the  problem  at 
hand,  the  so-called  EDMC  schemes.  The  conservation  laws  of  linear  and  angular  momenta, 
defined  in  this  discrete  setting 

B'node.  ^node 

ln+i  =  MAB  Vn+i  alld  3n+i  =  Xn+i  X  ,  (H-3.7) 

A,B= 1  A,B=  1 

(i.e.  ln+i  =  In  and  3n+ 1  =  3n  f°r  fext  =  0),  follow  easily  from  the  considered  mid-point 
approximation  (II. 3. 4),  as  long  as  the  B-bar  operator  satisfies  the  relations 

(II. 3. 8) 

for  a  constant  vector  c  e  Rn<lil"  (here,  (•)  denotes  the  global  vector  of  nodal  values  given  by 
(•)).  We  observe  that  the  conditions  (II. 3. 8)  are  the  discrete  counterparts  of  the  relations 
(II. 2. 3)  for  the  continuum  problem. 

The  group  motion  associated  to  the  relative  equilibria  of  the  discrete  equations  (II. 3.3)- 
(II. 3. 4)  were  obtained  in  Armero  &  Romero  [2001a]  and  are  given  by 

xAn  =  AntpheA  +  un  and  vAn  =  An  [f2e  x  ipheA  +  Ve\  ,  (II.3.9) 

for  fixed  vectors  f2e  and  Ve,  and  a  sequence  of  rotations  { An }  and  displacements  {iin } 
defined  recursively  by  the  relations 

An- |_i  =  An  CAY  At  SPIN  [J2e]  and  un+ i  =  un  +  At  An+ 1  Ve  ,  (II. 3. 10) 
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At  / 

1  +  —spin  [£2e] 


At 

I  -  —SPIN  [f2e 


1  -1 


for  An+ 1  :=  (A„  +  An+i)  /2  (not  a  rotation  in  general)  and  the  Cayley  transform 
CAY  At  SPIN  [fte 

for  the  discrete  relative  equilibrium  configuration  <//'  given  by  the  equation 


£  SO(flc lim)  i 

(II. 3. 11) 


M(f2ex(neX(p*  +  vS)+  j  hBTeS{<phe)dV  =  0,  (II. 3. 12) 

as  long  as  the  condition 


B*e  —  BeA  ,  i 


n+i  ’ 


(II. 3. 13) 


is  satisfied  along  the  group  motion  (II. 3. 9).  We  observe  that  equation  (II. 3. 12)  is  the 
discrete  counterpart  of  the  equilibrium  equation  (II. 2. 5).  Thus  the  algorithm  preserves  the 
relative  equilibria  of  the  system  as  long  as  the  B-bar  operator  satisfies  condition  (II. 3. 13), 
the  counterpart  of  relation  (II. 2. 6). 

Finally,  the  counterpart  of  the  energy  conservation/dissipation  equation  (II. 2. 7)  is 
obtained  as 


and 


Hn-\- 1  Hn 


1 


'B 


A V  dV  for  AV  =  S  :  -AC  -  AW  ,  (II.3.14) 

A 


H-n+i  r.'Pn+i  4  AIvn-\-i  T  /  dV  i  0  or  1 

2  JBh 


(II. 3. 15) 


for  the  discrete  system  (II.3.3)-(II.3.4),  as  long  as  we  have  the  relation 


B*  (dn- |_i  dn)  —  ^  AC 


(II. 3. 16) 


for  the  B-bar  operator  By  and  the  increment  of  the  assumed  strain  AC/2. 


Clearly,  the  interest  here  lies  in  the  discrete  dissipation  (II. 3. 14)  reproducing  exactly 
the  dissipation  of  the  continuum  system.  For  the  elastoplastic  model  (II.2.10)-(II.2.13), 
this  can  be  accomplished  by  considering  the  elastoplastic  decomposition  of  the  assumed 
deformation  gradient  F  =  FeFp  (see  Section  II. 4.1  below)  and  the  discrete  equations 


\  ([*”’]  A i  A C[F*]^j  -  AC')  =  A7  N+(S„q.)  , 
skew  f[C* J„+.  (F>+ 1  -  F£)  [F^j]  =  A7  Mwr(S„q.)  . 

0*  ■■=  0(5*,  q*)  <  0  ,  Ay  >  0  ,  Ay0*  =  0 


(11. 3. 17) 

(11. 3. 18) 

(11. 3. 19) 
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as  proposed  in  Armero  [2006],  with  q *  =  —ATi/Aa  and  the  stresses  given  by  = 
Fp  x  S*FpT  x  with 

n+  2  n+2 


S* 


+  2 


W(C‘n+1)  -  W(Cn)  -  2jg:  (|C']„+J)  :  AC' 


\cxl  1  AC'  [C'r;  ,(11.3.20) 


-1 


[c 


e'“!i  ACe  :  ACe 

o 


Jn+ 


Here  we  have  introduced  the  notation  [(-)]re+i  =  ((-)n.  +  (-)n+i)  /2.  The  formula  (II. 3. 20) 

corresponds  to  a  conserving  approximation  of  the  gradient  formula  (II. 2. 9)  such  that  S  : 
|ACe  =  AIFe.  It  is  a  modification  of  the  the  original  conserving  formula  proposed  in 
Gonzalez  [2000]  by  including  the  elastic  metric  [Ce]~+i ,  as  it  will  become  crucial  in  the 
construction  of  the  assumed  B-bar  operator  in  the  following  section.  Again,  the  exact 
plastic  dissipation  (including  exact  energy  conservation  for  an  elastic  step)  is  obtained 
by  the  return  mapping  algorithm  (II.3.17)-(II.3.19)  and  the  stress  formula  (II. 3. 20).  This 
situation  adds  the  desired  numerical  stability  to  the  algorithms  as  illustrated  with  the 
numerical  simulation  presented  in  Section  II. 5. 


A  variation  of  the  return  mapping  algorithm  (II.3. 17)- (II.3. 19)  that  also  imposes  ex¬ 
actly  the  isochoric  plastic  response  in  isochoric  plastic  models,  like  the  classical  J2-flow 
theory  of  metals,  can  be  found  in  Armero  &  Zambrana  [2007].  Similarly,  we  refer  to 
Armero  &  Romero  [2001a],  Armero  &  Romero  [2001b]  for  variations  of  the  stress  for¬ 
mula  (II. 3. 20)  that  incorporate  a  controllable  high-frequency  numerical  energy  dissipation 
to  handle  the  usual  high  numerical  stiffness  in  the  systems  of  interest. 


II. 4.  Conserving  assumed  strain  finite  element  methods 


II.4.1.  The  assumed  deformation  gradient  and  its  variations 


The  interest  here  is  the  development  of  assumed  strain  finite  element  methods  for 
the  locking-free  approximation  of  nearly  incompressible  material  models,  like  the  plas¬ 
ticity  models  outlined  above  combined  with  a  Mises-type  deviatoric  yield  surface,  while 
exhibiting  the  conservation/dissipation  laws  obtained  in  the  previous  section.  This  can 
be  accomplished  with  the  now  standard  scaled  deformation  gradient  (see  e.g.  Armero 
[2000],  SlMO  ET  AL  [1988]) 


F  n-\-i  — 


e 


n-\-i 


■J, 


n-\-i 


GRAD  [<£n+j]  for  Jn+i  :=  det  GRAD 


(11*4.1) 


(for  i  =  0,  or  1)  and  the  assumed  Jacobian  O  =  det  F  defined  by  the  weighted  average  at 
the  element  level 


en+i(x)  :=  r  T(X)  H1  [  r(Y)jn+i(Y)  dv , 


(II. 4. 2) 
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for  a  set  of  n$  element  interpolation  functions  T  (X)  and  the  corresponding  matrix 

H  :=  f  r(Y)  r  T(Y)  dV  E  ,  (II.4.3) 

over  a  generic  element  B]'; .  Typical  choices  are  no  =  1  and  T  =  1  in  combination  with  a 
bilinear  quadrilateral  or  trilinear  brick  element  (Q1/A0),  and  no  —  3  with  r  =  [1  £  rj]T  for 
the  isoparametric  coordinates  (£,  rf)  in  plane  problems  (and  similarly  for  3D)  in  combination 
with  quadratic  quadrilateral  or  triangular  elements  (Q2/A1  or  P2/A1). 

—  — T  — 

The  consideration  of  the  assumed  right  Cauchy-Green  tensor  Cn+i  =  Fn+iFn+i  leads 
to  the  variations  8Cn+i/2  =  Bn+i8dn+i  for 


-Dn+* 


e 


n-\-i 


J, 


n+i 


Bn+i  +  -  C(cp^+l)  ®  (g 


A 

n+i 


(II. 4. 4) 


defining  the  classical  B-bar  operator  Bn+i  in  terms  of  the  standard  (displacement)  lin¬ 
earized  strain  operator 


(vhn+i)T2  J V| 

(Vhn+i)l  N$  +  (rt+,)T2  Ni 


for  A  =  1,  nnode  , 


(II.4.5) 


for  plane  problems  (similarly  for  3D),  and  the  spatial  gradients 


9n+i  ■=  K?iGA  and  9n+i  -=rT  H1  f  r  gA+i  dV ,  (II.4.6) 

JBh 

for  the  material  gradients  GA  :=  GRAD[iV-4]  =  [iVj1  Ar^]T  (A  =  l,nnode )  as  used  in 
(II. 4. 5).  Here  we  consider  i  =  0,  1/2  or  1,  with  i  =  1/2  corresponding  to  the  evaluation  of 
the  different  quantities  above  in  the  mid-point  configuration  <^+1  =  +  (p!^+1  )/2. 

The  standard  choice  Bn+i  in  the  governing  equation  (II. 3. 3)  does  not  lead,  however, 
to  a  conserving  approximation.  The  conditions  (II. 3. 8)  can  be  easily  seen  to  be  satisfied, 
but  the  conditions  (II. 3. 13)  and  (II. 3. 16)  for  the  conservation  of  the  relative  equilibria  and 
energy  are  not.  This  situation  is  to  be  traced  to  the  spatial  gradients  gA  L  in  (II. 4. 6).  First, 
we  observe  that  during  the  group  motion  (II. 3. 9)  we  have  F„  =  AnFe  and  Fn+ 1  =  An+1  Fe 
for  the  equilibrium  deformation  gradient  Fe,  but 


I  —T  jp—Tf^iA  _  \—T  „A 

=  9? 


An+\  9e 


(II. 4. 7) 
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as  required  for  (II. 3. 13),  since  An+x  =  (An  +  An+i)  /2  is  not  a  rotation  in  general.  Simi¬ 
larly,  the  condition  (II. 3. 16)  is  not  satisfied,  in  part  because 


node. 


Y  •  idn+ 1  -  dn)  ± 


(Jn+1  Jn) 


A=1 


[J] 


n+h 


(II. 4. 8) 


for  [J]n_ |_i  =  (Jn  +  Jn+ 1)/2,  as  it  would  be  expected  from  its  continuum  counterpart 
( J  =  J  Vv  :  1  for  the  spatial  velocity  gradient  Vv). 


II. 4. 2.  A  new  conserving  B-bar  operator 

Faced  with  the  difficulties  observed  in  the  previous  section,  we  introduce  the  new 
modified  spatial  gradients 


9 


n+i 


rA1  r  - 1  \c\nU :  AC 

■■=  K+l  IC]'U  GA  +  2-lFy - ^  F_ 


[cidj  :  ac  Mdj 


n+i 


[crnUAc[c\~UG 


n+i 


for  A  =  1 ,  rinode  ,  and  their  assumed  counterparts 


_  1 

Sn+h  ~  J+i 


[0] 


rrH 


n+b 


'  LA  Wn+i9*+idV 


(II. 4. 9) 


(II. 4. 10) 


where  we  note  the  use  of  the  average  Jacobian  [J]n+i  =  (Jn  +  Jn+ 1 ) /2.  We  observe  that 
these  modified  spatial  gradients  do  satisfy,  by  construction,  the  relations 


n nod 


/  -<  c'n+2 
A= 1 


and 


ffn+i  ■  K+l  -  dn  )  = 


Jn+1  Jn 

“PI 


J^node 


and  Y  $ ™+h'  (dY  ~  d™) 


n+i 


A= 1 


—/I 


9n+ 1  =  An+ige  and  gn+ 1  =  An+ige  , 


@n+ 1 

(11. 4. 11) 

(11. 4. 12) 


along  the  group  motion  (II. 3. 9),  all  for  A  =  1  ,nnode-  We  can  observe  the  similarities  with 
the  stress  formula  (II. 3. 20)  and,  in  particular,  the  use  of  the  reference  (convected)  metric 
[C]n_ |_i  to  arrive  at  the  proper  transformation  properties  for  the  modified  spatial  gradients. 

With  these  definitions  at  hand,  we  introduce  the  following  new  B-bar  operator 


(II. 4. 13) 
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for  A  =  1,  nnode  where 

(  [Ml  ,  (('fo+i)  "((•)»)  fQr  M  J.  (.) 

£,(«)._  J  n+^  (■)„+!-(•)„  1OT  ^n+lT{)n, 

{  a  ([(-)]n+l)  for  (')„+!  =  (•)„  , 


(II. 4. 14) 


for  a  generic  quantity  (•)  and  exponent  a.  By  the  way,  we  note  that  the  formulas  (II. 3. 20) 
and  (II. 4. 9)  are  well-defined,  with  the  quotients  vanishing  when  C%+ ,  =  Ccn  and  = 

Cn,  respectively.  No  singularity  occurs. 

The  different  terms  in  the  expression  (II. 4. 13)  can  be  seen  to  be  second-order  approx¬ 
imation  of  the  variations  of  the  assumed  C.  Some  long  algebraic  manipulations  show  that 
this  new  B-bar  operator  satisfies  the  desired  conditions  (II. 3. 8),  (II. 3. 13)  and  (II. 3. 16).  In 
particular,  the  relation  (II. 3. 13)  for  the  relative  equilibria  is  satisfied  for  the  assumed  B-bar 
operator 


S"  +  3C' 


(b 


for  A  =  l,nnode  , 


(II. 4. 15) 


at  the  equilibrium  configuration  Hence,  the  new  B-bar  operator  (II. 4. 13)  leads  to  a 
fully  energy-momentum  assumed  strain  finite  element  formulation  when  combined  with 
the  EDMC  time-stepping  algorithms  considered  in  Section  II. 3,  obtaining  in  this  way  a 
(energy)  stable  formulation  that  avoids  volumetric  locking. 


II. 5.  Representative  numerical  simulations 

We  present  in  this  section  several  numerical  simulations  to  illustrate  and  evaluate 
the  performance  of  the  assumed  finite  elements  developed  in  this  paper.  In  particular, 
we  present  in  Section  II. 5.1  convergence  tests  showing  the  locking- free  response  of  the 
proposed  elements  in  the  incompressible  limit.  Section  II. 5. 2  evaluates  the  energy  and 
momentum  conservation  properties  in  time  for  a  plane  elastic  solid  in  free-fliglit,  including 
the  conservation  of  the  relative  equilibria.  Finally,  Section  II. 5. 3  evaluates  the  energy 
dissipation  and  momentum  conservation  properties  in  time  for  an  elastoplastic  solid  in  the 
general  three-dimensional  case. 


II. 5.1.  Cook’s  membrane  problem:  evaluation  of  the  locking-free  properties 

We  consider  in  this  section  the  classical  benchmark  problem  of  the  so-called  Cook’s 
membrane  for  the  evaluation  of  the  locking  properties  of  finite  elements  in  the  incompress¬ 
ible  limit.  Figure  II. 5.1  depicts  a  sketch  of  the  problem,  consisting  of  a  tapered  slab  fixed 
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FIGURE  II. 5.1  Cook’s  membrane  problem.  Problem  definition:  geometry,  boundary  condi¬ 
tions  and  loading.  Distances  in  mm.  Plane  strain  conditions  are  assumed  for  a  unit  thickness 
in  the  perpendicular  direction. 


TABLE  II. 5.1  Cook’s  membrane.  Material  parameters. 


Young  modulus 

E 

210 

GPa 

Poisson  ratio 

V 

0.2-0.4999 

Reference  density 

Po 

1.0T03 

kg/m3 

at  one  end  and  with  an  imposed  transversal  load  F(t)  at  the  opposite  end.  The  load  is 
applied  proportionally  in  time 

rn  =  P  r .  (n.5.i) 

accounting  for  the  dynamic  effects  as  it  is  the  interest  in  this  work.  It  is  a  dead  load, 
uniformly  distributed  along  the  reference  configuration  of  the  slab’s  right  edge.  The  solid 
is  assumed  to  be  at  rest  in  the  undeformed  configuration  at  t  =  0.  Plane  strain  conditions 
are  assumed. 

We  consider  the  compressible  Neo-Hookean  hyperelastic  model,  defined  by  the  poten¬ 
tial 

W{C )  =  -  (In  J)2  +  |  (trace  [C\  -  3)  -  [i  In  J  ,  (II.5.2) 

for  the  Lame  constants  2//,  =  E/(  1  +  u)  and  A  —  Ifw j (1  —  2v)  in  terms  of  reference 
the  Young  modulus  E  and  Poisson  ratio  v.  Table  II. 5.1  includes  the  values  of  these 
parameters  assumed  in  the  simulations  presented  here,  including  also  the  reference  density 
pa.  To  evaluate  the  performance  of  the  elements  in  the  incompressible  limit,  the  values 
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of  the  Poisson  ratio  v  =  0.2  and  0.4999  are  considered,  the  latter  leading  effectively  to  a 
quasi-incompressible  response. 

We  consider  the  three  characteristic  plane  elements:  the  two  quadrilateral  elements 
Q1/A0  and  Q2/A1  consisting  of  bi-linear  and  bi-quadratic  displacement  interpolations  with 
piece-wise  constant  and  linear  interpolations  for  the  volumetric  strain  interpolation,  re¬ 
spectively.  To  illustrate  the  development  of  a  triangular  element,  we  consider  the  P2+/A1 
element  with  7-node  quadratic  displacements  and  a  piece-wise  linear  interpolation  for  the 
volume.  The  simulations  presented  here  for  the  quadratic  elements  Q2/A1  and  P2+/A1  do 
consider  the  presence  of  an  internal  node  and  its  contribution  to  the  displacement,  velocity 
and  acceleration  interpolations  (and  hence,  with  mass  contributions).  Other  implemen¬ 
tations  where  these  contributions  are  treated  in  the  context  of  enhanced  strain  elements 
as  internal  bubbles  affecting  only  the  strains  can  be  easily  devised.  We  consider  a  consis¬ 
tent  mass  matrix  in  all  cases.  The  Q1/A0  and  Q2/A1  quads  consider  a  2  x  2  and  3x3 
(full)  Gauss  quadrature,  respectively,  whereas  the  quadratic  P2+/A1  element  considers 
the  standard  6  point  quadrature  on  triangles;  see  e.g.  Hughes  [1987]. 

The  interest  in  this  section  is  the  evaluation  of  the  locking  properties  of  the  spatial 
interpolations  defined  by  the  assumed  strain  elements  developed  here.  To  this  end,  we 
consider  simulations  for  different  levels  of  spatial  mesh  refinement,  all  with  a  fixed  time 
step.  We  run  the  problem  with  100  equal  time  steps  to  the  final  time  tf  =  10  ms  and  final 
load  F  =  100  k,N ,  and  measure  the  top  corner  vertical  displacement.  We  consider  regular 
structured  spatial  meshes  with  equal  number  of  nodes  per  side.  Figure  II. 5. 2  depicts  the 
solutions  obtained  with  the  three  different  elements  and  the  mesh  with  17  nodes  per  side 
for  quasi-incompressible  case  v  =  0.4999.  We  have  included  the  distribution  of  the  Mises 
(Cauchy)  stress  (i.e.  ||dev[crn+1]  ||  for  the  Cauchy  stress  rrn+]  :=  Fn+iS*F n+1/On+\ 
where,  recall,  On+i  =  det  )  superposed  to  the  deformed  spatial  configuration.  A 

good  agreement  can  be  observed  between  the  different  solutions.  We  note  the  presence 
of  large  displacements  and  strains.  The  energy-momentum  conserving  formula  (II. 3. 20)  is 
considered  for  the  evaluation  of  the  stresses  in  the  time-stepping  scheme  assumed  for  the 
temporal  integration  of  the  equation,  thus  leading  to  a  fully  conserving  approximation  in 
time.  We  evaluate  these  properties  in  the  following  sections. 

Figure  II. 5. 3  depicts  the  computed  top  corner  displacement  versus  the  number  of  nodes 
per  side  for  the  two  Poisson’s  ratios  of  v  =  0.2  and  v  =  0.4999.  A  fairly  good  performance 
is  observed  for  the  former  case  by  all  the  elements.  This  performance  deteriorates  drasti¬ 
cally  in  the  quasi-incompressible  case  v  =  0.4999  for  the  basic  displacement  finite  elements 
Ql,  Q2  and  P2+.  These  elements  simply  consist  of  the  evaluation  of  the  linearized  strain 
operator  Bn+ 1  in  (II. 4. 5)  at  the  mid-point  configuration  tn+ 1.  This  choice  leads  also  to 
the  sought  energy-momentum  conservation  properties  in  time  but  clearly  to  volumetric 
locking  in  the  incompressible  limit.  As  expected  the  performance  of  the  simple  bilinear 
displacement  Ql  quad  is  particularly  bad.  The  Q2/A1  and  P2+/A1  also  show  a  signifi¬ 
cant  improvement  over  their  displacement  counterparts  Q2  and  P2+  elements,  as  a  faster 
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FIGURE  II. 5. 2  Cook’s  membrane  problem.  Distribution  of  the  Mises  stress  (in  GPa )  super¬ 
imposed  to  the  deformed  configurations  at  the  final  time  t  =  10s  for  the  quasi-incompressible 
case  ( v  =  0.4999)  and  for  three  different  assumed  strain  elements:  the  bilinear  Ql/AO  and 
quadratic  Q2/A1  quads  and  the  quadratic  assumed  P2+/A1  triangle.  All  correspond  to  the 
meshes  with  17  nodes  per  side  and  follow  the  same  scale  shown. 
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v  =0.4999 


FIGURE  II. 5. 3  Cook’s  membrane  problem.  Convergence  plots  for  two  different  Poisson’s 
ratios  ( v  =  0.2  and  v  =  0.4999)  and  six  different  finite  elements.  The  removal  of  the  locking 
in  the  quasi-incompressible  case  v  =  0.4999,  to  be  contrasted  with  the  corresponding  basic 
displacement  element,  is  to  be  noted. 
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FIGURE  II. 5. 4  Two-dimensional  elastic  solid  in  free  flight.  Problem  definition:  geometry  and 
initial  conditions.  Distances  in  m.  Plane  strain  conditions  are  assumed  for  a  unit  thickness  in 
the  perpendicular  direction. 


TABLE  II. 5. 2  Two-dimensional  elastic  solid  in  free-flight.  Material 
parameters  for  the  elastic  arms. 


Young  modulus 

E 

40 

GPa 

Poisson  ratio 

V 

0.45 

Reference  density 

Po 

8.6-103 

kg/ vp? 

convergence  is  observed  for  limit  values  obtained  with  finer  meshes.  When  comparing  the 
plots  presented  in  Figure  II. 5. 3,  we  conclude  that  the  proposed  assumed  treatment  of  the 
volumetric  strain  leads  to  the  desired  locking-free  response  in  volume. 

II. 5. 2.  Two-dimensional  solid  in  free  flight:  evaluation  of  the  conservation 

properties  in  time  for  the  elastic  case,  including  the  relative  equilibria 

We  evaluate  in  this  section  the  conservation  properties  in  time  of  the  newly  developed 
B-bar  formulation  for  nonlinear  elastic  problems.  We  still  consider  plane  strain  conditions, 
leaving  the  consideration  of  three-dimensional  problems  for  the  next  section.  The  interest 
here  is  the  confirmation  of  the  conservation  of  the  energy  and  the  preservation  of  the  mo¬ 
mentum  conservation  laws,  and  the  associated  relative  equilibria,  as  shown  in  Proposition 
1.3.1. 

To  this  purpose,  we  consider  the  plane  solid  depicted  in  Figure  II. 5. 4  consisting  of  a 
rigid  cylindrical  core  and  two  elastic  arms.  This  figure  also  shows  the  finite  element  mesh 
considered  in  all  the  numerical  simulations  presented  here.  It  consists  of  72  quadrilat¬ 
eral  assumed  strain  Q1/A0  elements,  defined  by  bilinear  displacements  and  the  piece- wise 
constant  volumetric  strain.  The  consistent  mass  matrix  (II. 3. 5)  is  considered. 

The  elastic  arms  follow  the  compressible  Neo-Hookean  response  defined  by  the  stored 
energy  function  (II. 5. 2).  Table  II. 5. 2  includes  the  values  of  the  material  parameters  as¬ 
sumed  in  the  simulations.  In  particular,  we  consider  a  Poisson  ratio  of  v  =  0.45  for  the 
problem  to  be  close  to  the  incompressible  limit,  as  evaluated  in  the  previous  section.  The 
rigid  core  has  been  modeled  by  simply  increasing  the  Young  modulus  to  100  times  that  of 
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the  arms. 

The  solid  is  given  the  initial  velocity  defined  by  the  nodal  values 

vA(0)  =  f2  e3  x  xA(0)  A=  1,  nnode  ,  (II.5.3) 

for  an  initial  angular  velocity  f2,  the  unit  vector  e3  corresponding  to  the  thickness  direction 
(i.e.  the  direction  perpendicular  to  the  plane  of  the  problem)  and  the  initial  nodal  positions 
xA  (0)  =  XA  of  the  reference  mesh  in  Figure  II. 5. 4  in  a  Cartesian  reference  system  with 
origin  at  the  center  of  the  rigid  core.  The  solid  is  left  then  to  evolve  freely  in  the  plane 
{ei,  62},  with  e\  defined  by  the  axis  of  the  arms  in  their  initial  configuration.  Figure  II. 5. 5 
shows  the  motion  computed  with  the  newly  proposed  B-bar  method  (II. 4. 13)  for  an  initial 
angular  velocity  of  Q  =  0.230  rad/ms  and  a  constant  time  step  of  At  =  0.5  ms.  We 
observe  that  the  rigid-body  velocity  distribution  (II. 5. 3)  is  lost  as  the  solid  deforms  from 
this  initial  velocity  and  undeformed  configuration.  We  have  included  the  distribution  of 
the  Mises  Cauchy  stress  in  the  arms,  defined  as  indicated  in  the  previous  section,  over  the 
deformed  configuration  of  the  solid.  Large  strains  and  displacements  (rotations)  around 
the  center  (by  symmetry)  can  be  observed. 

Given  the  assumed  free- flight  conditions,  the  linear  momentum  lh  and  the  angular 
momentum  jh  should  be  conserved  along  the  motion.  In  particular,  jh  =  jh  e3  for  a 
scalar  jh  in  this  plane  setting  and  lh  =  0  for  the  assumed  initial  conditions.  Similarly, 
the  total  energy  Hh  in  (II. 3. 14)  should  also  be  conserved  for  the  assumed  elastic  response. 
These  conservation  properties  are  confirmed  for  the  new  B-bar  method  (II. 4. 13)  as  shown 
in  Figure  II. 5. 6,  depicting  the  exact  conservation  of  the  angular  momentum  jh  and  energy 
Hh  along  the  motion.  The  different  components  of  the  linear  momentum,  not  shown, 
are  also  exactly  conserved.  The  kinetic  and  potential  energies,  see  equation  (II. 2. 7),  are 
also  included  in  the  energy  plot  with  their  sum,  the  total  energy,  showing  again  its  exact 
conservation.  These  results  confirm  Proposition  1.3.1. 

To  illustrate  the  need  of  the  new  B-bar  operator  for  the  exact  conservation  of  the 
energy  (that  is,  to  satisfy  the  relation  exactly)  we  consider  the  standard  (non-conserving) 
B-bar  operator  (II. 4. 4)  evaluated  at  the  mid-point  Bn+i.  As  discussed  above,  this  choice 
does  conserve  the  linear  and  angular  momenta  but  not  the  energy.  This  situation  is 
confirmed  by  the  numerical  simulations.  Figure  II. 5. 7  shows  the  evolution  of  the  total 
energy  for  this  case  and  compares  it  with  the  new  conserving  B-bar  method  (II. 4. 13),  for 
two  different  time  steps  (At  =  0.5  and  0.3  ms).  The  conserving  stress  formula  (II. 3. 20) 
is  maintained,  so  the  only  difference  is  the  B-bar  operator.  The  lack  of  the  conservation 
of  the  energy  for  the  standard  B-bar  operator  (II. 4. 4)  is  apparent,  even  though  again  a 
conserving  stress  formula  is  used. 

It  is  interesting  to  note  that  the  solutions  in  Figure  II. 5. 7  begin  to  show  high-frequency 
content  in  the  long  term,  leading  to  a  non-convergence  of  the  simulation  for  the  non¬ 
conserving  B-bar  operator  at  a  higher  energy  content.  It  must  be  noted  that  these  failures 
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FIGURE  II. 5. 5  Two-dimensional  elastic  solid  in  free  flight:  solution  for  an  angular  velocity 
17  =  0.230  rad/ms  (At  =  0.5  ms).  The  distribution  of  the  Mises  stress  (in  GPa)  in  the  arms  is 
shown  superimposed  to  the  solid’s  deformed  configuration  at  different  times. 
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FIGURE  II. 5. 6  Two-dimensional  elastic  solid  in  free  flight:  solution  for  an  angular  velocity 
17  =  0.230  rad/ms  (A t  —  0.5  ms).  Evolution  in  time  of  the  angular  momentum  (left)  and 
the  energy,  including  the  kinetic,  potential  and  total  energies  (right),  computed  with  the  new 
conserving  assumed  strain  finite  element  in  combination  with  an  energy-momentum  conserving 
time-stepping  scheme.  The  exact  conservation  of  the  total  energy  and  the  angular  momentum 
is  to  be  noted. 
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FIGURE  II. 5. 7  Two-dimensional  elastic  solid  in  free  flight.  Solution  for  an  angular  velocity 
17  =  0.230  rad/ms.  Evolution  in  time  of  the  total  energy  computed  with  the  new  conserving 
B-bar  operator  (II.4.13)  and  the  basic  mid-point  B-bar  operator  (II. 4. 4),  for  different  time 
steps. 
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FIGURE  II. 5. 8  Two-dimensional  elastic  solid  in  free  flight:  group  motion  associated  with 
the  relative  equilibrium  at  an  angular  velocity  ft  =  0.230  rad,/ ms  (At  =  0.5  ms).  The 
distribution  of  the  Mises  stress  (in  GPa)  in  the  arms  is  shown  superimposed  to  the  solid’s 
deformed  configuration  at  different  times.  An  exact  rotation  can  be  observed  for  all  times. 


Non-equilibrium  solution 


Relative  equilibrium 


FIGURE  II. 5. 9  Two-dimensional  elastic  solid  in  free  flight.  Evolution  in  time  of  the  radius 
for  nodes  A,  B  and  C  (see  Figure  II. 5. 4)  for  the  solutions  starting  from  the  undeformed  con¬ 
figuration  (left)  and  the  relative  equilibrium  (right),  both  for  the  same  initial  angular  velocity 
ft  =  0.230  rad/ms  (At  =  0.5  ms). 
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of  convergence  also  affect  fully  energy  conserving  formulations  (although  they  are  observed 
to  fail  to  converge  at  a  later  time),  being  a  direct  consequence  of  the  high  numerical  stiffness 
of  the  problem.  We  refer  to  Armero  &  Romero  [2001a],  Armero  &  Romero  [2001b] 
for  a  complete  discussion  of  this  issue  as  well  as  the  formulation  of  EDMC  time-stepping 
algorithms  that  exhibit  the  controllable  high-frequency  energy  dissipation  required  to  han¬ 
dle  these  situations.  Briefly,  this  is  accomplished  with  a  modified  stress  formula  (II. 3. 20), 
so  the  new  stress  approximation  S'*  satisfies  the  energy  relation 

S*  :  ^  AC  =  AW  +  Vnum  ,  (II. 5. 4) 

for  the  numerical  dissipation  Vnum  >  0,  affecting  the  high-frequency  response  without 
perturbing  the  momentum  conservation  properties,  being  rigorously  non-negative  and  fully 
controllable  (i.e.  T>num  =  0  is  an  option,  recovering  the  original  fully  energy  conserving 
approximation  as  a  particular  case).  The  important  fact  for  the  discussion  here  is  that 
the  B-bar  operator  still  needs  to  satisfy  the  relation  (II.3.16)  for  the  numerical  dissipation 
Vnurn  in  (II. 5. 4)  to  appear  in  the  global  energy  evolution  equation,  and  hence  for  the 
scheme  to  be  rigorously  energy  dissipative  (conservative  as  a  particular  case).  Therefore, 
the  consideration  of  these  EDMC  schemes  requires  also  the  use  of  the  new  conserving  B- 
bar  operator  developed  in  this  paper  with  the  basic  operator  (II. 4. 4)  still  prone  to  exhibit 
numerical  instabilities.  The  original  references  Armero  &  Romero  [2001a],  Armero  & 
Romero  [2001b]  considered  only  displacement  based  finite  element  formulations,  which 
do  satisfy  these  conditions  as  discussed  above  but  lead  to  locking  elements. 

To  complete  the  confirmation  of  Proposition  1.3.1,  we  study  the  conservation  of  rela¬ 
tive  equilibria  for  the  problem  at  hand.  To  this  purpose,  we  solve  the  relative  equilibrium 
equation  (II. 3. 12),  obtaining  the  equilibrium  deformation  pf/(X)  for  a  given  angular  ve¬ 
locity  Qe.  Since  this  deformation  is  defined  up  to  a  translation  and  a  rotation,  and  given 
the  symmetry  in  the  motions  considered  here  (with  the  center  of  mass  remaining  fixed  at 
the  origin),  the  boundary  conditions  imposed  while  solving  this  equation  consist  of  fixing 
the  node  at  this  center  and  the  transversal  displacement  of  Node  C  (see  Figure  II. 5. 4)  at 
the  tip  of  an  arm  (transversal  to  the  axis  of  the  arm),  hence  fixing  the  free  rigid  rotation. 

We  then  compute  the  dynamic  solution  with  the  initial  conditions  defined  by  the 
computed  equilibrium  deformation  (i.e.  xA (0)  =  pl/(XA))  and  the  velocity  distribution 
(II. 5. 3)  for  the  equilibrium  angular  velocity  Qe  based  on  this  equilibrium  configuration. 
If  the  relative  equilibria  are  preserved,  the  group  motion  given  by  the  solution  (II. 3. 9), 
corresponding  to  a  rigid  rotation  in  the  deformed  configuration  ip^(X)  with  a  constant 
angular  velocity  Qe  (note  that  ve  =  0,  i.e.,  no  translation),  must  be  recovered. 

The  new  B-bar  operator  does  obtain  this  solution,  hence  confirming  the  analyses 
leading  to  Proposition  1.3.1.  Figure  II. 5. 8  depicts  the  solution  obtained  in  this  case  with 
an  angular  velocity  i?e  =  0.230  rad/ ms  (At  =  0.5  ms).  The  deformation  at  t  =  0  ms  shows 
the  computed  equilibrium  configuration  cp1/  (X)  with  the  distribution  of  the  Mises  (Cauchy) 
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stress  in  the  arms;  see  Section  II.5.1.  The  same  configuration,  with  the  very  same  stress 
distribution,  is  obtained  at  later  times  but  rotated.  The  purely  rigid  rotational  character 
of  this  motion  is  also  confirmed  by  the  plots  in  Figure  II. 5. 9,  showing  the  evolution  in 
time  of  the  radial  coordinate  from  the  center  of  the  rigid  core  of  Nodes  A,  B  and  C  along 
the  arms;  see  Figure  II. 5. 4.  These  radii  remain  constant  at  all  times  at  the  initial  value 
corresponding  the  initial  deformed  configuration  c p^(X).  Note  that  these  radii  do  not 
correspond  to  the  values  in  the  initial  undeformed  configuration,  except  for  Node  A  since 
this  node  is  on  the  rigid  central  core. 

To  illustrate  better  this  equilibrium  solution,  we  have  included  in  Figure  II. 5. 9  the 
radii  of  the  same  nodes  for  the  dynamic  solution  depicted  in  Figure  II. 5. 5.  In  contrast  with 
the  equilibrium  solution,  this  non-equilibrium  solution  starts  from  the  undeformed  configu¬ 
ration  shown  in  Figure  II. 5. 4.  The  same  initial  angular  velocity  of  Qe  =  0.230  rad/ms  has 
been  considered  in  both  cases.  The  variation  of  these  radii  illustrate  the  non-equilibrium 
character  of  this  solution  in  contrast  with  the  solution  depicted  in  Figure  II. 5. 8.  Note  again 
that  Node  A  is  on  the  rigid  central  core.  We  also  note  that  these  two  solutions,  although 
defined  by  the  same  initial  angular  velocity,  correspond  to  different  angular  momentum 
and  energy,  both  conserved  in  the  equilibrium  and  non-equilibrium  solutions  but  different 
given  the  initial  conditions  for  each. 

II. 5. 3.  Three-dimensional  solid  in  free  flight:  evaluation  of  the  conservation/ 
dissipation  properties  in  time  for  the  elastoplastic  case 

The  purpose  of  this  last  numerical  example  is  to  evaluate  the  conservation  properties 
in  time  of  the  formulation  developed  in  this  paper  in  the  context  of  elastoplastic  problems 
and  in  the  three-dimensional  case.  As  noted  above,  elastoplastic  problems  are  of  particular 
interest  since  the  isochoric  plastic  flow  in  typical  applications  involving  plastic  models  like 
1/2 -flow  theory  leads  effectively  to  quasi-incompressible  problems,  thus  requiring  special 
locking-free  finite  element  interpolations  as  developed  in  this  paper.  We  also  consider  a 
three-dimensional  problem  in  contrast  with  the  plane  problem  of  the  previous  section.  In 
this  way,  we  consider  the  assumed  strain  brick  Q1/A0  consisting  of  a  trilinear  interpolation 
of  the  displacement  with  a  piece-wise  constant  volume  over  the  element.  The  standard 
2x2x2  Gauss  quadrature  is  considered  in  the  evaluation  of  all  the  element  arrays. 

The  solid  under  consideration  consists  of  a  stiff  heavy  inner  ring  with  three  equally 
spaced  flexible  arms  in  the  configuration  depicted  in  Figure  II. 5. 10.  Details  of  the  geomet¬ 
ric  dimensions  are  shown  in  this  figure.  The  flexible  arms  follow  a  finite  strain  elastoplastic 
model  of  ./2-flow  theory  based  on  a  multiplicative  decomposition  of  the  deformation  gradi¬ 
ent  F  =  FeFp.  Briefly,  we  consider  the  plastic  evolution  equations  (II.2. 10)- (II.2. 13)  for 
the  von  Mises  yield  function 


0  =  1 1  dev  [t] 


^  y(a)  <  0  , 


(II. 5. 5) 
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FIGURE  II. 5. 10  Three-dimensional  solid  in  free  flight.  Problem  definition:  geometry  and 
initial  conditions.  Distances  in  m. 


for  the  Kirchhoff  stress  r  :=  FeSFe  ,  its  (spatial)  deviatoric  part  dev  [r]  and  the  corre¬ 
sponding  Euclidean  tensor  norm  (i.e.  ||A||2  :=  A  :  A),  with  the  plastic  flow  vectors 


JV*= 

os 


i,C* 


dev  [r]  reT 

~F  lidSTT 


n<j,  = 


d(j) 

dq 


SC « 


and  M\yp  =  0 


(IL5.6) 

(that  is,  an  associated  plastic  model  with  no  plastic  spin) ,  after  noting  that  we  have 


1 1  dev  [r]  ||2  =  CeDEVc-i 


DEVc-i 


Ce 


for 


DEVc->  [s]  :=  S  -  i  (s  :  C*)  C'“  =  F'dev  [r]  F‘ 


(II.5.7) 


(IL5.8) 


thus  confirming  the  assumed  functional  dependence  for  the  yield  function  (II. 5. 5)  and  flow 
vectors  (II. 5. 6).  The  linear  isotropic  hardening  law 


y(a)  =  ay  +  H  a  ,  (II. 5. 9) 

and  corresponding  hardening  potential  7i(a)  =  H  a2/ 2  (so  q  =  — dTL/da  =  —Ha  = 
<jy  —y(a)).  for  the  initial  uniaxial  yield  limit  ay  and  linear  hardening  modulus  H.  has  been 
assumed. 

The  elastic  part  (II. 2. 9)  is  defined  by  the  Hencky’s  hyperelastic  potential 

We(Ce)  =  ^  (In  Je)2  +  “||  lnCe||2  ,  (II. 5. 10) 

Z  t: 

T 

for  the  elastic  right  Cauchy-Green  tensor  Ce  =  Fe  Fe,  the  elastic  Jacobian  Je  =  det Fe, 
and  the  Lame  constants  A  and  //,  as  in  the  previous  section.  Table  II. 5. 3  includes  the 
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TABLE  II. 5. 3  Three  dimensional  solid  in  free-flight.  Material  param¬ 
eters  for  the  elastoplastic  outer  arms. 


Young  modulus 

E 

10 

GPa 

Poisson  ratio 

V 

0.2 

Initial  uniaxial  yield  limit 

Gy 

975 

GPa 

Linear  hardening  modulus 

H 

200 

MPa 

Reference  density 

Po 

2.0-103 

kg/m3 

material  parameters  assumed  for  the  arms.  The  inner  ring  is  stiff  and  heavy,  following  the 
elastic  Hencky’s  model  (II. 5. 10)  with  a  reference  Young  modulus  100  Earrns  and  reference 
density  10  pQ  of  the  corresponding  values  for  the  arms,  and  with  the  same  Poisson  ratio  of 
v  =  0.2. 

We  refer  to  Armero  [2006],  Armero  &  Zambrana  [2007]  and  references  therein  for 
complete  details  on  the  elastoplastic  model  and  the  development  of  EDMC  return  mapping 
algorithms  for  their  integration.  In  particular,  we  consider  the  EDMC  return  mapping  algo¬ 
rithm  developed  in  Armero  [2006] ,  and  outlined  above,  in  combination  with  the  invariant 
stress  formula  (II. 3. 20)  in  terms  of  the  elastic  Cauchy-Green  tensor  Ce  for  the  evaluation 
of  the  stresses  (II. 2. 9)  in  its  elastic  part.  The  resulting  global  time-stepping  algorithms 
show  the  exact  plastic  dissipation  (including  exact  energy  conservation  in  elastic  steps) 
and  the  exact  conservation  laws  of  linear  and  angular  momenta.  The  successful  extension 
of  these  conservation/dissipation  properties  to  the  assumed  strain  elements  needed  for  fi¬ 
nite  element  simulations  of  problems  considering  this  type  of  plastic  model  requires  the 
developments  presented  in  the  current  paper.  We  also  note  that,  even  though  we  present 
results  here  for  the  particular  model  given  by  (II.5.5)-(II.5.10),  the  algorithms  considered 
here  are  completely  general,  with  formulas  like  (II. 3. 20)  and  (II.3.17)-(II.3.19),  applying  to 
general  isotropic  or  anisotropic  models,  associated  or  even  non-associated  plastic  models. 
As  noted  in  above,  it  is  precisely  the  flexibility  in  handling  different  material  models  that 
motivates  the  development  of  the  assumed  strain  framework  considered  here. 

The  simulations  are  started  by  considering  the  initial  velocity  given  by  the  rigid-body 
distribution  (II. 5. 3),  in  terms  of  the  angular  velocity  17  around  the  axis  of  symmetry  e-$ 
shown  in  Figure  II. 5. 10,  for  the  inner  ring  only.  No  loading  or  boundary  conditions  are 
imposed,  thus  leading  to  the  free-flight  of  the  solid  in  later  times.  In  these  conditions, 
the  three  components  of  both  the  linear  and  angular  momenta  are  conserved  as  it  is  the 
total  energy  for  elastic  cases.  A  positive  energy  dissipation  is  to  be  observed  if  the  plastic 
response  of  the  arms  is  activated. 

To  illustrate  the  fully  conserving  properties  in  the  elastic  case  in  the  considered  three- 
dimensional  setting,  we  first  consider  an  initial  angular  velocity  of  17  =  0.100  rad/ ms  , 
which  leads  to  a  fully  elastic  response.  The  computed  solution  for  a  constant  time  step  of 
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FIGURE  II. 5. 11  Three-dimensional  solid  in  free  flight:  solution  for  an  angular  velocity 
fi  =  0.100  rad/ms  leading  to  a  purely  elastic  response  of  the  solid.  The  distribution  of  the 
Mises  (Cauchy)  stress  (in  GPa )  in  the  arms  is  shown  superimposed  to  the  solid’s  deformed 
configuration  at  different  times. 
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FIGURE  II. 5. 12  Three-dimensional  solid  in  free  flight:  solution  for  an  angular  velocity 
17  =  0.100  rad /ms  leading  to  a  purely  elastic  response  of  the  solid.  Evolution  in  time  of 
the  angular  momentum  components  (left)  and  total  energy  (right)  computed  with  the  new 
conserving  assumed  strain  finite  element  in  combination  with  an  energy-momentum  conserving 
time-stepping  scheme.  The  exact  conservation  of  the  total  energy  and  the  three  components  of 
the  angular  momentum  is  to  be  noted. 


At  =  1.0  ms  is  depicted  in  Figure  II. 5.11.  The  distribution  of  the  Mises  stress  in  the  arms 
is  again  shown  over  the  deformed  configuration  of  the  whole  solid  for  different  times.  No 
activation  of  the  plastic  response  is  observed  for  these  initial  conditions,  thus  expecting 
a  full  conservation  of  the  total  energy.  This  conservation  property  is  confirmed  by  the 
evolution  plots  shown  in  Figure  II. 5. 12.  The  potential  and  kinetic  energies  are  shown  as 
they  vary  in  time  with  their  sum  (the  total  energy)  being  exactly  constant  during  the 
simulation.  Figure  II. 5. 12  shows  also  the  evolution  of  the  three  components  of  the  angular 
momentum  jh,  confirming  again  its  exact  conservation  at  all  times.  Note  that,  for  the 
considered  initial  condition,  the  only  non-zero  component  is  the  component  j%.  The  same 
situation  applies  to  the  linear  momentum,  not  shown,  with  the  three  components  vanishing 
given  the  assumed  initial  conditions  (lh  =  0). 

The  increase  of  the  initial  angular  velocity  to  17  =  0.200  rad/ms  leads  to  the  solution 
depicted  in  Figure  II. 5. 13.  This  initial  velocity  results  in  the  development  of  plastic  strains 
in  the  arms.  Figure  II. 5. 13  shows  the  distribution  of  the  equivalent  plastic  strain  a  at 
different  early  stages  of  the  motion.  The  large  deformation  and  strains  are  to  be  noted. 
A  constant  time  step  of  At  =  1.0  ms  has  also  been  considered.  Figure  II. 5. 14  depicts  the 
evolution  of  the  total  energy  and  the  axial  (non-zero)  component  of  the  angular  momen¬ 
tum  for  this  case.  The  exact  conservation  of  the  momentum  is  again  obtained,  with  the 
total  energy  showing  a  strictly  positive  energy  dissipation  during  plastic  steps  and  exact 
conservation  during  elastic  ones. 

These  results  confirm  the  exact  conservation  properties  of  the  assumed  strain  inter- 
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FIGURE  II. 5. 13  Three-dimensional  solid  in  free  flight:  early  stages  of  the  solution  for  an 
angular  velocity  f2  =  0.200  rad/ ms  leading  to  an  elastoplastic  response  of  the  solid.  The 
distribution  of  the  equivalent  plastic  strain  is  shown  superimposed  to  the  solid’s  deformed 
configuration  at  different  times. 
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FIGURE  II. 5. 14  Three-dimensional  solid  in  free  flight:  solution  for  an  angular  velocity  17  = 
0.200  rad/ms  leading  to  an  elastoplastic  response  of  the  solid.  Evolution  in  time  of  the 
angular  momentum  jQ  (left)  and  total  energy  (right)  for  the  new  conserving  assumed  strain 
finite  element  in  combination  with  an  EDMC  time-stepping  scheme  and  the  standard  (non¬ 
conserving)  assumed  strain  element  with  the  classical  trapezoidal  rule.  The  latter  exhibits  an 
instability  leading  to  the  termination  of  the  computation  with  a  non-physical  growth  of  the 
energy. 


polations  developed  in  this  work  in  combination  of  EDMC  time-stepping  algorithms  for 
the  integration  of  the  equations  in  time.  This  situation  is  to  be  contrasted  with  standard 
schemes  currently  available  in  the  literature.  As  an  illustration  we  have  included  in  Figure 
II. 5. 14  the  solution  obtained  with  the  classical  trapezoidal  rule  (i.e.  Newmark  j3  =  1/4, 
7  =  1/2);  see  e.g.  SlMO  [1998]  for  an  implementation  in  the  finite  strain  plastic  case  of 
interest.  In  particular  the  simulation  considers  the  exponential  return  mappings  that  can 
be  found  developed  in  this  last  reference.  The  spatial  discretization  considers  the  classical 
(non-conserving)  B-bar  approximation  of  Section  (sub:NCBbar)  evaluated  at  the  end  of  the 
time  step  tn+ 1,  as  required  by  this  scheme.  The  lack  of  conservation/dissipation  properties 
for  both  the  angular  momentum  and  energy  become  apparent  in  the  solutions  depicted  in 
Figure  II. 5. 14.  In  particular,  we  observe  the  nonlinear  unstable  character  of  the  scheme  as 
the  unbounded  growth  of  energy  forces  to  stop  the  simulation  at  some  time.  The  presence 
of  this  unphysical  energy  growth,  even  in  the  context  of  a  dissipative  elastoplastic  as  con¬ 
sidered  here,  is  to  be  noted  and  contrasted  with  the  exact  energy  conservation/dissipation 
properties  of  the  methods  developed  in  this  work. 


II. 6.  Concluding  remarks 

We  have  presented  in  this  appendix  a  new  framework  for  the  development  of  assumed 
strain  finite  elements  for  nearly  incompressible  problems  in  the  finite  deformation  dynamic 
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range  that  preserves  the  conservation/dissipation  properties  in  time  of  energy  and  mo¬ 
menta  characteristic  of  recently  developed  energy-momentum  and  EDMC  time-stepping 
algorithms.  The  proposed  formulation  relies  on  the  proper  definition  of  the  linearized  strain 
operator  (or  B-bar  operator)  accounting  not  only  for  the  spatial  interpolations  but  also 
for  the  temporal  approximation  of  the  governing  equations.  Standard  existing  treatments, 
with  the  B-bar  operator  defined  through  the  standard  variation  of  the  assumed  strains, 
have  been  shown  to  destroy  the  energy  conservation,  thus  leading  to  difficulties  in  the  form 
of  instabilities  with  an  energy  growth.  After  identifying  the  general  conditions  that  this 
operator  must  satisfy  to  lead  to  conserving  approximations  of  the  energy,  the  linear  and 
angular  momenta  and  the  associated  relative  equilibria,  we  have  undertaken  the  formula¬ 
tion  of  actual  finite  elements  for  plane  and  three-dimensional  problems.  The  conservation 
properties  in  time  have  been  proven  rigorously,  and  evaluated  numerically  with  several 
representative  simulations  in  nonlinear  elastodynamics  and  dynamic  finite  strain  plastic¬ 
ity,  illustrating  also  the  locking-free  character  of  the  final  elements  in  the  incompressible 
limit  of  interest. 

A  main  advantage  of  the  considered  assumed  strain  framework  is  its  flexibility  in  terms 
of  the  involved  constitutive  models.  Indeed,  the  kinematic  considerations  defining  the  new 
B-bar  operator  are  completely  independent  of  the  considered  material  model,  this  being 
elastic  or  inelastic,  isotropic  or  anisotropic,  with  coupled  or  uncoupled  volumetric  and 
deviatoric  parts.  This  situation  is  to  be  contrasted  with  the  more  involved  implementation 
of  mixed  formulations  for  general  models;  see  e.g.  SiMO  et  al  [1988].  We  have  only 
addressed  in  this  paper  the  case  of  volumetric  locking,  without  considering  shear  locking. 
The  latter  is  not  as  a  severe  problem  for  the  type  of  problems  considered  in  this  paper, 
namely,  continuum  problems;  see  e.g.  Armero  [2000],  Armero  [2004].  However,  the 
assumed  strain  methods  developed  here  for  the  volume  response  can  be  combined  with 
enhanced  strain  formulations  to  address  shear  locking,  as  developed  in  Armero  [2000]  for 
static  problems.  We  are  currently  exploring  these  extensions  for  the  dynamic  problems  of 
interest  here. 


References 


Armero,  F.  [2000]  “On  the  locking  and  stability  of  finite  elements  in  finite  deformation 
plane  strain  problems,”  Computers  &  Structures,  75:261-290. 

Armero  F  [2004]  “Assumed  strain  finite  element  methods,”  in  Finite  Element  Methods: 
1970’s  and  beyond,  ed.  by  LP  Franca,  TE  Tezduyar  and  A  Masud,  CIMNE,  Barcelona. 


F.  Armero 


72 


ARMERO,  F.  [2006]  “Energy-dissipative  momentum-conserving  time-stepping  algorithms 
for  finite  strain  multiplicative  plasticity,”  Comp.  Meth.  Appl.  Mech.  Eng.,  195:4862- 
4889. 

Armero,  F.  &  Romero,  I.  [2001a]  “On  the  formulation  of  high-frequency  dissipative 
time-stepping  algorithms  for  nonlinear  dynamics.  Part  I:  Low  order  methods  for  two 
model  problems  and  nonlinear  elastodynamics,”  Comp.  Meth.  Appl.  Mech.  Eng., 
190:2603-2649. 

Armero,  F.  &  Romero,  I.  [2001b]  “On  the  formulation  of  high-frequency  dissipative 
time-stepping  algorithms  for  nonlinear  dynamics.  Part  II:  High  order  methods,”  Comp. 
Meth.  Appl.  Mech.  Eng.,  190:6783-6824. 

ARMERO,  F.  &  Romero,  I.  [2003]  “Energy-dissipative  momentum-conserving  time-stepping 
algorithms  for  the  dynamics  of  nonlinear  Cosserat  rods,”  Comp.  Mech.,  31:3-26. 

Armero,  F.  &  Zambrana,  C.  [2007]  “Volume-preserving  energy-dissipative  momentum- 
conserving  algorithms  for  isochoric  multiplicative  plasticity,”  Comp.  Meth.  Appl.  Mech. 
Eng.,  196:4130-4159. 

Crisfield,  M.  &  Shi,  J.  [1994]  “A  co-rotational  element /time- integration  strategy  for 
non-linear  dynamics,”  Int.  J.  Num.  Meth.  Eng.,  37:1897-1913. 

Gonzalez,  O.  [2000]  “Exact  energy- momentum  conserving  algorithms  for  general  models 
in  nonlinear  elasticity,”  Comp.  Meth.  Appl.  Mech.  Eng.,  190:1763-1783. 

Hughes,  T.J.R.  [1980]  “Generalization  of  selective  integration  procedures  to  anisotropic 
and  nonlinear  media,”  Int.  J.  Num.  Meth.  Eng.,  15:1413-1418. 

Hughes  TJR  [1987]  The  Finite  Element  Method,,  Prentice-Hall. 

Marsden,  J.E.  [1992]  Lectures  on  Mechanics,  London  Mathematical  Society  Lecture  Note 
Series,  174:  Cambridge  University  Press. 

MENG,  X.N.  &  LAURSEN,  T.A.  [2002]  “Energy  consistent  algorithms  for  dynamic  finite 
deformation  plasticity,”  Comp.  Meth.  Appl.  Mech.  Eng.,  191:1639-1675. 

Nagtegaal,  J.C.;  Parks,  D.M  &  Rice,  J.R.  [1974]  “On  numerically  accurate  finite 
element  solutions  in  the  fully  plastic  range,”  Comp.  Meth.  Appl.  Mech.  Eng.,  4:153- 
177. 

Noels,  L.;  Stainier,  L.  &  Ponthot,  J.P.  [2004]  “An  Energy-Momentum  Conserving 
Algorithm  for  Non-Linear  Hypoelastic  Constitutive  Models,”  Int.  J.  Num.  Meth.  Eng., 
59:83-114. 

Romero,  I.  &  Armero,  F.  [2002]  “Numerical  integration  of  the  stiff  dynamics  of  ge¬ 
ometrically  exact  shells:  an  energy-dissipative  momentum- conserving  scheme,”  Int.  J. 
Num.  Meth.  Eng.,  54:1043-1086. 


Final  Report.  FA9550-05-1-0117 


73 


SlMO  JC  [1998]  “Numerical  Analysis  of  Classical  Plasticity,”  Handbook  for  Numerical 
Analysis,  Volume  VI,  ed.  by  P.G.  Ciarlet  and  J.J.  Lions,  Elsevier,  Amsterdam. 

SlMO,  J.C.  &  TARNOW,  N.  [1992]  “The  discrete  energy-momentum  method,  conserving 
algorithms  for  nonlinear  elastodynamics,”  ZAMP,  43:757-793. 

SlMO,  J.C.;  Posbergh,  T.A.  &  Marsden,  J.E.  [1991]  “Stability  of  relative  equilibria. 
Part  II:  Application  to  nonlinear  elasticity,”  Arch.  Rational  Mech.  Anal.,  115:61-100. 

SlMO,  J.C.;  Taylor,  R.L.  &  Pister,  K.  [1988]  “Variational  and  projection  methods 
for  the  volume  constraint  in  finite  deformation  elastoplasticity,”  Comp.  Meth.  Appl. 
Mech.  Eng.,  57:177-208. 


